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ABSTRACT 


The  objectives  of  this  program  were  the  study  of  the  chemistry  of  YBCO  thin  film  deposition  by 
MOCVD,  and  the  model-based  control  of  MOCVD  reactors.  Model  reduction  techniques  were 
applied  to  chemical  vapor  deposition  (CVD)  of  YBCO  thin  films.  This  w^k  lias  paralleled  some 
of  the  work  of  Goodwin  et  al  at  Caltech  sponsored  under  the  DARPA  VIP  Phase  I  program,  but 
with  significantly  different  approach  and  emphasis. 

Under  this  program,  the  above  chemistry  was  investigated  from  first  principles  using  quantum 
chemistry  computations  (led  by  MIT).  The  kinetic  and  thermodyn^ic  information  obtained 
from  these  studies  were  used  to  generate  a  kinetic  mechanism,  which  was  then  coupled  to 
transport  models  for  fluid  flow,  heat  and  mass  transfer  in  CVD  reactors.  These  coupled  reaction- 
transport  models  were  used  to  derive  reduced  order  models  for  process  control.  In  addition,  ef 
was  directed  toward  model  reduction  (led  by  SC  Solutions  and  Princeton  University). 

SC  Solutions  implemented  the  CVD  modeling  techniques  on  a  commercially  available  Cro/CW 
software.  Because  of  its  wide  acceptance  in  the  semiconductor  industry,  we  selected 
from  CFDRC  as  the  prototype  platform.  In  particular,  SC  Solutions  devel(^ed  a  model  o  a 
system  used  by  Superconductor  Technologies  Inc.  (STI)  for  deposition  of  YBCO  high  thin  fi  ms 

using  CFD-ACE. 
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Executive  Summary 

sc  Solutions,  MIT,  and  Princeton  have  developed  computer  simulation  tools  and  control 
strategies  for  production  of  yttrium-barium-copper  oxide  (YBaiCusOT-x  or  YBCO)  thin  films  with 
high  temperature  superconducting  (HTS)  properties,  using  metal-organic  chemical  vapor 
deposition  (MOCVD)  process.  HTS  thin  films  have  applications  in  high  speed-low  power 
computing  and  signal  processing  electronics,  high  performance  microwave  device^  magnetic 
sensors,  and  in  optoelectronic  systems  for  communications  and  signal  analysis.  MOCVD  otters 
the  potential  for  growth  under  highly  oxidizing  conditions,  for  large  area  deposition,  and  high 
throughput.  Our  research  on  CVD  modeling  involves  concurrent  development  of  first  principles 
chemical  mechanisms  with  reactor  scale  modeling  to  provide  a  more  realistic  description  of 
MOCVD  reactors.  Earlier,  very  little  was  known  of  the  chemical  mechanism  involved  in 
decomposition  of  metal  precursors  and  oxide  formation  both  in  gas-phase  and  at  the  wafer 
surface.  There  had  not  been  any  systematic  investigation  of  control  issues  for  better  control  of 
growth  rates,  oxide  stoichiometry,  and  film  thickness  uniformity.  Our  teani  has  successfully 
addressed  the  above  issues  through  fundamental  quantum-mechanical  calculations,  reactor-scale 
physical  process  model,  model  order  reduction,  and  developing  model-based  control  strategies. 
These  simulation  tools  and  the  results  obtained  from  the  studies  provides  a  clear  understanding  of 
the  chemical  mechanism,  species  transport,  and  film  growth.  This  understanding  enables  design 
and  implementation  of  optimized  controllers  and  that  meet  both  process  specifications  as  well  as 
run-to-run  repeatability  which  is  essential  for  large-scale  production.  The  following  summanzes 
the  research  results  in  more  detail. 

First,  principle  quantum  chemistry  calculations  have  been  performed  to  probe  the  reaction 
pathways  in  the  decomposition  of  metal  P-diketonates,  used  as  precursors  for  this  CVD  process. 
The  metal  p-diketonates  can  be  represented  by  the  general  formula  M(RiC(0)CHC(0)R2)n  where 
M  is  the  metal  (M  =  Cu,  Y),  Ri  and  R2  are  the  substituents  on  the  ligand  (Ri=R2=CH3  for  acac 
and  t-butyl  for  dpm)  and  n  denotes  the  number  of  ligands  (n=2  for  Cu,  n=3  for  Y).  There  is  ve^ 
little  information  available  on  these  compounds  which  makes  the  choice  of  the  method  difficult. 
Different  methods  and  basis  sets  have  been  tested  and  evaluated  by  comparing  model  predictions 
against  experimental  data  for  small  copper  compounds  having  similar  bonds  as  the  copper 
p-diketonates.  Based  on  our  validation  studies,  hybrid  density  functional  theory  method, 
particularly  the  implementation  using  the  three  parameter  Becke  exchange  and  the  Lee-Yang-Pair 
correlation  functionals  (B3LYP),  was  chosen.  For  copper,  we  used  the  effective  core  potentials 
by  Hay  and  Wadt  [1-3]  with  a  Neon  core  with  the  associated  double  zeta  basis  set.  The 
Dunning/Huzinga  full  double  zeta  (D95)  basis  set  was  used  for  all  other  elements. 

Mean  bond  dissociation  energies  (BDE)  of  the  M-0  bond  (M  =  Cu,  Y)  and  the  scission  of  the 
substituents  Ri,  R2  (CH3  and  t-butyl)  from  the  ligand  have  been  determined.  Relative  magnitudes 
of  the  bond  strengths  suggest  the  break-up  of  the  chelate  ring  as  the  initial  step,  followed  by  the 
cleavage  of  the  t-butyl  (or  CH3)  group  from  the  precursor.  This  is  accompanied  by  further  break¬ 
up  of  the  entire  organic  structure.  The  Cu-0  bond  rupture  seems  to  be  the  rate-controlling  step. 
Based  on  these  observations,  a  kinetic  mechanism  has  been  formulated.  Rate  parameters  have 
been  estimated  for  the  rate-controlling  step  and  they  correlate  well  with  the  rate  parameters 
obtained  from  the  overall  decomposition  process  using  mass-spectrometry. 

Structures  obtained  for  Y(acac)3  are  in  good  agreement  with  experimental  data.  The  Y-0  bond 
strength  is  higher  compared  to  the  Cu-0  bond  strength  in  Cu(acac)2.  Morokuma’s  “our  own  n- 
layered  integrated  molecular  orbital  and  molecular  mechanics”  (ONIOM)  method  [4]  allows  the 
treatment  of  different  parts  of  a  molecule  at  different  levels  of  theory,  and  has  been  used  to  study 
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the  precursor,  Y(dpm)3.  The  bulky  t-butyl  t-butyl  groups  have  been  treated  using  molecular 
mechanics  and  the  central  rings  are  treated  using  DFT  methods.  Y(acac)3  has  been  used  as  a 
model  compound  for  Y(dpm)3  in  some  cases.  The  decomposition  is  initiated  by  the  cleavage  of 
the  t-butyl  group  followed  by  ligand  pruning  to  give  Y(dpm)2.  These  results  are  in  agreement 
with  experimental  observations. 


HTS  MOCVD  models  were  developed  for  the  2D  version  of  the  commercial  STl  reactor  using 
CFD-ACE,  and  deposition  rates  and  film  thickness  uniformity  of  YBCO  films  were  obtained. 
The  deposition  rates  compare  favorably  with  the  published  literature  and  the  experimental  results 
from  STl.  Results  were  obtained  on  Principal  Orthogonal  Decomposition  (POD)-based 
compression  of  two  kinetic  schemes.  A  general  control  structure  was  developed  for  MOCVD 
reactor  control  systems.  An  innovative  run-to-run  controller  scheme  was  developed  that  enables 
efficient  convergence  to  the  desired  stoichiometry  region  for  MOCVD  systems  using  a  multi-step 
process  control  method. 
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Accomplishments 

•  Developed  a  systematic  methodology  and  formulated  strategy  to  understand  the  chemical 
mechanism  involved  in  MOCVD  of  YBCO  thin-films. 

•  Undertook  fundamental  investigation  using  density  functional  theory  for  quantum- 
mechanical  calculations  to  determine  most  probable  pathways  for  precursor  decomposition 
and  oxide  formation  leading  to  YBCO  thin-film  growth.  These  calculations  have  greatly 
increased  understanding  of  precursor  decomposition.  They  studies  resolved  controversies  in 
the  literature,  namely  the  existence  of  some  gas-phase  reactions  as  part  of  the  decomposition 
mechanism. 

•  This  quantum  chemistry  methodology  allows  determination  of  rate  constants,  etc.,  in  absence 
of  difficult  and  costly  experiments,  and  allows  determination  of  the  effects  of  impurities 

•  The  above  approach  for  elucidating  the  chemical  mechanism  involved  in  the  decomposition 
of  such  large  molecules  is  a  pioneering  effort! 

•  A  two-dimensional  MOCVD  reactor-scale  model,  developed  using  an  approximate  chemical 
mechanism,  has  been  used  for  calculating  growth  rate,  deposition  uniformity,  and  oxide 
stoichiometry  at  wafer  surface.  Results  comparable  to  those  in  the  literature  and  those 
recorded  at  an  industrial  reactor  (STI,  Santa  Barbara,  CA).  The  model  was  refined  using  rate 
constants  from  the  DFT  calculations. 

•  Control  architecture  has  been  developed  for  both  run-to-run  and  real-time  feedback. 
Response  surface  maps  were  used  to  develop  model-based  controllers  for  better  uniformity, 
and  for  run-to-run  variations  due  to  bubbler  temperature  drift,  etc.  The  developed  run-to-run 
control  scheme  allows  control  of  stoichiometries  of  Ba,  Y,  and  Cu  of  the  YBCO  thin  film. 
The  scheme  is  very  innovative  as  the  process  model  can  be  used  in  lieu  of  costly  experiments 
to  provide  quick  convergence  to  the  desired  stoichimetry  region. 
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1.  Introduction  and  Background  on  HTS  and  CVD  Modeling 


HTS 

High  temperature  superconductivity  (HTS)  has  applications  in  technologies  such  as  high-speed, 
low-power  computing  and  signal  processing  circuits,  high  performance  microwave  devices, 
magnetic  sensors,  and  optoelectronic  systems  for  communications  and  signal  analysis.  These 
materials  have  several  potential  applications,  including  magnetic  film  detection,  fast  computing 
devices,  interconnects  for  multichip  modules,  bolometric  detectors,  and  microwave  devices. 
Complex  oxide  materials  are  used  in  thin  film  form  to  make  such  devices  and  circuits. 

CVD  Modeling 

The  focus  of  this  research  effort  was  on  chemical  mechanism  generation  and  reduction  for 
metalorganic  chemical  vapor  deposition  (MOCVD)  of  high  temperature  superconducting  (HTS) 
films.  As  a  technique  for  fabrication  of  HTS  films,  MOCVD  offers  the  potential  for  growth  under 
highly  oxidizing  conditions,  for  large  area  deposition,  and  high  throughput.  However,  the  process 
is  highly  complex  involving  flow,  heat  transfer,  and  gas-phase  and  surface  chemical  kinetics  at 
high  temperatures.  The  decomposition  reactions  of  the  precursors  are  crucial  to  the  process. 
Typical  precursors  for  the  growth  of  YBa2Cu307.x  films  are  p-diketonates,  e.g.,  Ba(dpm)2 
(dpm=3D  dipivaloylmethanate)  or  Ba(hfa)2=B7tet  (hfa=3D  hexafluoroacetylacetonate,  tet=3D 
tetraglyme).  These  precursors  are  typically  difficult  to  transport  and  their  decomposition 
mechanisms  are  poorly  characterized.  The  current  understanding  of  MOCVD  of  HTS  films  is 
limited  to  conditions  leading  to  mass  transport  or  kinetically  limited  growth  regimes.  Semi- 
empirical  models  based  on  assumptions  of  a  few  gas-phase  reaction  steps  followed  by  overall 
surface  reactions  have  been  proposed,  but  the  chemical  mechanisms  underlying  MOCVD  is  not 
understood.  An  improved  understanding  of  the  basic  chemistry  could  lead  to  process 
improvements  and  would  provide  the  basis  for  model-based  process  control  algorithms  through 
the  use  of  nonlinear  model  reduction  techniques.  These  issues  were  addressed  in  this  program. 
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2.  Gas-Phase  Chemistry  of  High  Temperature 
Superconductors 


Introduction 


Many  fundamental  challenges  in  the  growth  of  high  temperature  superconductor  (HTS) 
films  by  metalorganic  chemical  vapor  deposition  (MOCVD)  are  associated  with  the  properties  of 
the  precursor  compounds  [5-8].  Metal  p— diketonates  are  the  most  widely  used  precursors  and 
have  the  general  formula  M(RiC(0)CHC(0)R2)n,  where  M  is  the  metal  atom,  Ri  and  R2  are  the 
substituents  in  the  ligand  and  n  is  the  number  of  ligands  in  the  complex.  Some  of  the  common 
substituents  employed  in  HTS  film  deposition  are  Ri  =  R2  =  CH3  (acetylacetonate,  acac),  Ri  —  R2 
=  C(CH3)3  (dipivaloymethanate,  dpm)  and  Rj  =  R2  =  C(CF)3  (hexafluroracetylacetonate,  hfac). 
The  general  structure  of  a  copper  P-diketonate  is  shown  in  Figure  1.  MOCVD  of  HTS  films  is  a 
complex  process  involving  both  gas-phase  and  surface  reactions.  Deposition  of  HTS  films  of 
definite  composition  requires  the  knowledge  of  decomposition  products  and  kinetic  p^ameters. 
Choice  of  the  substituents  on  the  ligand  influences  the  volatility  and  thermal  stability  of  the 
precursor  molecule.  The  precursor  molecules  are  t3q)ically  solids  at  room  temperature.  They  are 
either  vaporized  before  transporting  to  the  deposition  chamber,  or  are  dissolved  in  an  organic 
solvent  and  then  delivered  via  aerosol  delivery  to  the  deposition  chamber.  The  latter  process  is 
complicated  by  combustion  of  the  organic  solvent  molecules.  The  reactant  gases  are  mixed  with 
an  oxidant  gas  such  as  O2,  N2O  [9]  or  O3  [10]  and  undergo  decomposition  and  oxidation  over  the 
heated  substrate  to  deposit  a  film  of  the  desired  material. 


Figure  1.  Schematic  of  the  Cu  p-diketonate. 


YBa2Cu307  is  the  most  widely  studied  HTS  material.  Common  precursors  in  the 
MOCVD  of  YBCuO  films  are  Y(dpm)3,  Ba(dpm)2  Cu(acac)2  and  Cu(dpm)2.  There  are  very 
limited  experimental  data  for  the  thermolysis  and  oxidation  of  the  metal  P-diketonate  compounds. 
Only  recently,  kinetic  studies  using  mass  spectroscopy  and  FTIR  spectroscopy  have  begun  to 
appear. 


The  gas-phase  reactions  in  the  CVD  chamber  can  have  a  marked  influence  on  film 
properties  such  as  uniformity,  impurity  incorporation  and  deposition  rates.  For  ex^ple,  films 
grown  using  the  hfac  ligand  are  found  to  have  a  high  fluorine  contamination.  Similarly,  using 
N2O  as  the  oxidant  gas  (in  place  of  O2)  lowers  the  substrate  temperature  requirements,  but  also 
decreases  the  growth  rate  by  a  factor  of  two  [11]. 
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We  have  used  quantum  chemistry  methods  to  explore  the  different  reaction  pathways  to 
gain  an  understanding  of  the  structure  and  reactivity  of  the  precursor  molecules.  The  attention 
particularly  focused  on  determination  of  the  initial  step  in  the  decomposition  process. 
Information  on  the  subsequent  reaction  steps  involving  the  organic  products  can  then  be  obmined 
from  the  extensive  literature  available  in  the  field  of  combustion  chemistry  [12],  The  kinetic  and 
thermodynamic  information  obtained  from  these  studies  can  be  used  to  generate  a  kinetic 
mechanism,  which  can  then  be  coupled  to  transport  models  for  fluid  flow ,  heat  and  mass  transfer 
in  CVD  reactors.  These  coupled  reaction-transport  models  can  be  used  to  derive  reduced  order 
models  for  process  control. 

Computational  Details: 

Gaussian98  suite  of  programs  [13]  was  used  for  all  calculations,  except  the  validation 
calculations  where  Gaussian94  and  MOLPRO  [14]  were  also  used.  The  hybrid  DPT  approach 
was  used  with  the  three  parameter  Becke  exchange  [15],  and  the  Lee- Yang-Parr  correlation 
functionals  [16]  (B3LYP).  For  copper,  we  used  the  effective  core  potentials  by  Hay  and  Wadt  [1- 
3]  with  a  Neon  core  with  the  associated  double  zeta  basis  set,  and  the  Dunning/Huzinga  full 
double  zeta  (D95)  basis  set  for  all  other  elements.  This  combination  is  referred  to  as  ECP II  here. 
The  choice  of  this  method/basis  set  combination  is  the  result  of  a  detailed  testing  of  different 
methods  and  basis  sets  and  comparing  the  results  with  experimental  data.  This  is  described  in  the 
subsequent  paragraphs. 

Validation  of  Methods  and  Basis  Sets 

The  application  of  quantum  chemistry  (QC)  methods  to  the  precursors  employed  in  CVD 
of  high  temperature  superconductors  is  a  challenging  task  due  to  the  large  number  of  electrons 
(many  of  them  core),  electron  correlation  and  increased  importance  of  relativistic  effects  for  large 
nuclei.  To  control  the  computational  time  needed  for  the  calculations,  two  parameters  can  be 
varied:  calculation  method  and  eventually  the  basis  set. 

Traditionally,  three  types  of  methods  have  been  distinguished:  empirical,  semi-empirical 
and  ab-initio  calculations.  For  molecules  of  the  size  used  in  the  current  project,  semi-empirical 
and  ab-initio  methods  are  the  most  common  ones.  In  case  of  the  ab-initio  methods,  two 
subgroups  can  be  divided,  i.e.  density  functional  theory  (DFT)  and  ‘traditional’  ab-initio  methods 
like  Hartree-Fock  (HF),  MoUer-Plesset  perturbation  theory  (MPx)  or  configuration  interaction 
(Cl).  Depending  on  the  level  of  accuracy,  the  ‘traditional’  methods  include  different  amounts  of 
the  correlation  energy  in  the  system,  however,  they  also  require  more  computational  time  with 
increasing  accuracy.  This  is  due  to  the  unfavorable  scaling  with  the  number  of  electrons  in  the 
system  for  the  modest-  to  high-level  methods  (like  Cl).  DFT  methods,  on  the  other  hand, 
generally  include  the  correlation  energy,  but  scale  much  more  favorable.  However,  the  crux  lies 
in  the  suitable  choice  of  the  density  functionals  applied,  which  is  still  not  fully  understood  due  to 
the  relatively  new  nature  of  the  DFT  method  itself. 

During  the  last  two  years,  hybrid  methods  like  Morokuma’s  ONIOM  [17,  18]  have  been 
developed,  which  can  couple  any  of  the  three  types  of  methods  in  a  single  calculation  on  a  single 
system.  In  this  case,  high  accuracy  methods  might  be  used  for  atoms  at  the  reaction  center  and 
computational  less  expensive  ones  for  atoms  in  the  periphery. 

The  common  type  of  a  basis  set  treats  every  single  electron  in  the  system,  which  is  to  be 
described.  Using  these  all-electron  basis  sets,  going  down  one  row  in  the  periodic,  the 
computational  cost  can  increase  by  up  to  a  factor  of  ten.  A  better  approach  is  the  use  of  the 
effective  core  potentials  (ECPs),  where  the  electronically  less  important  core  electrons  are 
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replaced  with  a  small  number  of  potentials,  resulting  in  substantial  savings  in  time,  memory  and 
disk  space  requirements.  As  a  further  advantage,  they  often  include  Darwin  and  mass-velocity 
relativistic  effects.  However,  the  decision  which  electrons  are  put  into  the  core  and  which  are 
explicitly  described  by  basis  functions  is  often  crucial  to  the  performance  of  the  ECPs.  For 
instance,  in  the  case  of  transition  metals  (electronic  configuration  [...]  (n-l)^"  (n-l)p®  (n-l)£/*  os^ 
np°)  the  correlation  between  (n-1)  s,p  electrons  and  the  (n-l)(/  electrons  can  be  so  strong  that  the 
whole  (n-1)  shell  has  to  be  excluded  from  the  core. 

Finding  a  combination  of  method  and  basis  sets,  which  does  not  exceed  a  reasonable 
amount  of  computer  resources,  but  still  give  sufficiently  accurate  results,  is  one  of  the  first  tasks 
in  any  computational  project.  In  the  current  project,  the  challenge  was  that  experimental  and 
computational  data  for  the  investigated  systems  are  sparse  and  are  not  sufficient  to  check  the 
accuracy  of  our  calculations.  Therefore,  an  evaluation  of  different  methods/basis  set  combination 
was  done  on  smaller  systems,  which  include  the  bond  types  involved  in  the  current  systems. 

Semiempirical  -  PM3 

The  large  size  of  the  precursor  molecules  would  suggest  the  use  of  semi-empirical 
models,  if  applicable.  PM3™  [19],  a  semi-empirical  method,  parameterized  for  transition  metals 
has  been  implemented  in  Spartan  [20]  and  its  use  was  initially  explored  for  bis-acetylacetonate 
copper  (H).  Its  structure  has  been  determined  to  have  copper  in  square-planar  coordination  and 
the  two  ligands  lying  in  one  plane  [21,  22].  However,  a  geometry  optimization  using  the  PM3™ 
method  does  not  agree  at  ^1  with  the  experimental  results,  but  rather  suggests  a  contorted 
structure  where  the  carbon  atoms  of  one  ligand  are  tipped  up  relative  to  the  rest  of  the  molecule 
by  almost  exacdy  90°.  Therefore,  this  approach  was  considered  inadequate  for  the  systems 
investigated. 

Ab  initio  -  DFT 

As  the  systems  to  be  investigated  in  this  project  are  relatively  large,  we  did  not  consider 
‘traditional’  ab  initio  methods  of  a  sufficient  accuracy  (MPx,  Cl)  to  be  reasonable  in  terms  of 
computational  time.  Instead,  we  were  focusing  on  the  DFT  methods,  checking  different 
functionals.  These  functionals  (and  combinations,  respectively)  were: 

(1)  Becke’s  correlation  functional  [15]  together  with  the  correlation  functional  of  Lee,  Yang  and 
Parr  [16]  (BLYP) 

(2)  Becke’s  3-parameter  hybrid  functional  [23]  together  with  the  correlation  functional  of  Lee, 
Yang  and  Parr  [16]  (B3LYP) 

(3)  Slater  exchange  fonctional  [24-26]  together  with  Vosko,  Wilk  and  Nusair’s  correlation 
functional  [27]  (LSDA) 

Configuration  interaction  (Cl)  and  Coupled  Cluster  methods  (e.g.  CCSD(T))  were  used  only  for 
comparison,  as  they  are  considered  to  treat  correlation  with  a  high  accuracy. 

For  the  basis  functions,  two  all-electron  basis  set  (6-31 IG  with/without  diffuse 
(+)/polarization  (*)  functions  and  Ahlrich’s  triple-^  valence  [tzv])  were  compared  with 
combinations  of  Hay  and  Wadt’s  effective  core  potentials  (Cu)  and  appropriate  basis  sets  for  all 
other  atoms  (ECP1/ST03-G  [ECPl]  and  ECP2/D95  [ECP2]).  In  some  cases,  the  Stuttgart  ECP 
has  also  been  used  [28].  There  had  been  doubts  about  the  applicability  of  Hay  and  Wadt’s  ECPs, 
along  with  their  associated  valence  basis  sets,  for  DFT  calculations.  In  a  study  on  CuCH2  [29], 
the  Hay-Wadt’ s  ECPs  have  been  checked  for  DFT  work  and  the  concerns  dismissed. 

The  systems  used  for  the  evaluation  of  the  different  bonds  were: 
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general  electronic 
structure 

Cu:  ionization  potential  [IP], 
excitation  energy  [Eexec] 

Cu-H 

Cu-H 

Cu-C 

Methyl-copper  (CuMe), 

Ethyl-copper  (CuEt) 

Cu-0 

Hydroxy-copper  (CuOH), 

Copper-dimethylate  (Cu(OMe)2), 

Copper  monoxide  (CuO) 

CC,  CO  and  CH  bonds  have  not  been  evaluated,  because  they  have  been  tested  extensively 
throughout  the  literature  and  found  to  be  described  with  good  accuracy  by  the  methods/basis  sets 
in  question.  Data  for  the  validation  studies  are  given  in  Appendix  I. 


Cu 

The  calculations  showed  that  the  excitation  energy  ->  can  only  be  computed 
reasonably,  if  the  correlation  energy  is  included  to  a  high  extent  (as  in  CCSD(T)).  DFT 
calculations  (like  B3LYP)  do  not  even  give  qualitatively  correct  results'.  Not  much  difference  has 
been  seen  for  the  different  ECPs. 

The  ionization  potential  is  much  less  sensitive  to  correlation  and  therefore  DFT  methods 
are  giving  reasonably  accurate  results.  B3LYP  was  found  to  perform  extremely  well,  with  errors 
not  more  than  6(2)  kcal/mol  for  all-electron  (ECPs),  except  if  6-3 IIG  is  used  without 
polarization/diffuse  functions. 

CuH 

For  CuH  a  short  evaluation  of  data  showed  the  influence  of  the  correlation  on  the 
dissociation  energy,  if  B3LYP  is  used.  In  this  case,  ECPl  did  not  perform  well,  due  to  the 
relatively  large  core.  ECP2  and  6-31 1+G**,  on  the  other  hand,  showed  about  the  same  quality  as 
highly  correlated  Cl  calculations  (including  core-core  correlation). 

CuMe,  CuEt 

For  methyl-copper  some  experimental  data  (vibrational  frequencies,  estimate  of 
dissociation  energy)  are  available,  which  makes  it  a  good  testing  candidate.  The  drawback, 
however,  is  the  pronounced  anharmonicity  of  the  methyl  group,  which  will  show  as  a 
disagreement  between  the  experimental  set  of  vibrational  frequencies  (including  anharmonicity) 
and  the  calculated  harmonic  ftindamentals  (without  anharmonicity,  as  the  name  implies). 

The  calculated  dissociation  energies  (exp.  60  kcal/mol)  for  CuMe  (6-31 1+G**,  ECPl, 
ECP2)  are  too  high  for  LSDA  calculations,  which  is  usual  for  this  method.  BLYP  calculations 
are  within  10  kc^mol  of  the  experimental  estimate,  with  the  exception  of  ECP2,  which  performs 
surprisingly  poorly  (35  kcal/mol  difference).  B3LYP  calculations  showed  good  agreement,  with 
6-31 1+G**  and  ECP2  about  the  same  value  (52.4  kcal/mol  and  53.7  kcal/mol).  ECPl  gives  a 
slightly  higher  value  (59.1  kcal).  With  CuEt  it  has  been  shown  that  the  diffuse  function  (+)  in  the 
6-3 IIG  basis  set  is  crucial  to  the  accuracy.  Polarization  functions  (**)  on  the  other  hand,  only 
slightly  change  the  calculated  energy  difference  (0.7  kcal/mol). 

To  further  evaluate  the  quality  of  the  basis  sets,  the  basis  superposition  error  was 
calculated  for  6-31 1+G**  (-5.6  kcal/mol),  ECPl  (8.0  kcal/mol)  and  ECP2  (-2.9  kcal/mol).  The 
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relative  high  value  for  the  ECPl  shows  that  this  ECP/basis  set  combination  is  not  flexible  enough 
to  correctly  describe  the  electron  distribution  in  the  system.  The  values  for  6-31 1+G**  and  ECP2 
are  considered  to  be  reasonable,  especially  because  relaxation  of  the  methyl  radical  after 
dissociation  has  not  been  taken  into  the  correction.  Correction  of  the  dissociation  energy  using 
the  BSSE  gives  6-31 1+G**  and  ECP2  values  within  3.5  kcal/mol  of  the  experimental  estimate 
and  extremely  close  to  the  values  of  Ziegler  [30]  and  Boehme  [31].  The  calculated  frequencies 
were  compared  with  the  experimental  data  set  for  6-31 1+G**,  ECPl  and  ECP2.  The  CH 
vibrations  are  usually  about  2-10%  overestimated  by  our  calculations  (6-31 1+G**  and  ECP2), 
but  come  out  12-20%  too  high,  if  ECPl  is  used.  For  the  Cu-Me  vibrations  (stretch  and  tilt)  50- 
60%  overestimation  seem  to  be  common,  with  the  Cu-Me  tilt  being  more  off.  Particularly  bad 
performing  are  B3LYP/ECP1  and  all  basis  sets  under  LSDA  (up  to  75%  overestimation).  These 
relatively  high  errors  are  probably  due  to  anharmonicity  as  well  as  method/basis  set  insufficiency. 
However,  no  information  can  be  obtained  on  how  much  each  of  these  factors  account  for. 

For  CuEt  similar  calculations  have  been  performed,  but  could  not  be  compared  to  any 
experimental  data.  A  factor  f  =  (o(B3LYP/X)/(o(B3LYP/6-311G)  has  been  introduced  to 
compare  the  calculational  data.  However,  this  factor  is  close  to  1  for  all  basis  sets  under  B3LYP 
and  almost  all  frequencies  (except  for  the  Cu  -  C=C  stretch).  No  particular  difference  regarding 
the  ECPl  could  be  seen. 

Cu-O  bond 

To  check  the  applicability  of  the  methods  regarding  the  most  important  bond  in  the 
present  systems,  different  models  could  be  thought  of.  The  one  being  chemically  similar  are 
copper  aUcoxides  like  Cu(OMe)2.  Unfortunately  these  compounds  are  unstable  towards 
hydrolysis  and  therefore  experimental  data  are  sparse.  Additionally,  the  experimental  data 
available  are  obtained  from  solid  state  where  copper  has  square-planar  coordination. 

These  differences  are  believed  to  lead  to  the  considerable  problems  faced  when  trying  to 
calculate  this  system  as  a  gas  phase  molecule.  For  B3LYP/6-311G**  and  B3LYP/ECP2 
geometry  optimizations  converged  and  harmonic  frequencies  were  compared  with  the 
experimental  ones.  Again,  an  overestimation  of  20  —  35  %  is  seen,  however,  data  did  not  allow 
for  a  further  analysis. 

Instead,  CuOH  and  CuO  were  used  to  check  for  the  calculation  of  the  Cu-O  bond, 
although  these  bonds  have  no  chemical  similarity  to  the  Cu-O  bond  in  the  investigated  systems. 
CuOH  geometries,  dissociation  energies  and  frequencies  were  compared  with  ECPl  and  ECP2 
using  B3LYP  and  SD-CI  methods.  It  again  becomes  clear  that  B3LYP/ECP2  is  close  in  accuracy 
to  all-electron  calculations  and  SD-CI  calculations,  but  ECPl  is  not  applicable. 

The  calculated  dissociation  energies  for  CuO  show  again  an  overestimation  by  at  least 
20-30  kcal/mol  for  LSDA  calculations  and  good  agreement  for  B3LYP/6-311+G**  and 
B3LYP/ECP2.  However,  in  contrast  to  the  calculations  at  CuMe/CuEt,  the  diffuse  function  does 
not  improve  the  accuracy  if  no  polarization  function  are  present. 

Summary  of  gas-phase  chemistry 

The  above  validation  of  basis  sets  and  methods  show  that  out  of  the  3  tested  functional 
combinations  for  DFT  calculations,  B3LYP  gives  generally  the  best  results  compared  with  the 
experimental  data.  Regarding  the  basis  set,  either  an  all-electron  triple-zeta  basis  set  (6-31 IG)  or 
a  Ne  core  ECP  (ECP2)  can  be  used  with  good  results.  However,  for  the  all-electron  basis  set,  it  is 
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important  to  include  both  diffuse  and  polarization  functions  as  they  influence  different  kinds  of 
bonds  in  different  ways. 

Structures 

Structural  parameters  obtained  from  quantum  chemistry  calculations  were  compared 
against  experimental  data  for  Cu(acac)2  [21,  22]  and  Y(dpm)3  [32],  obtained  from  electron 
diffraction  studies,  and  were  found  to  be  in  good  agreement.  Both  Cu(acac)2  and  Cu(dpm)2  have 
copper  in  a  square-planar  coordination  with  Cu-0  bond  length  of  ~  1.95  A.  Y(acac)3  has  a 
traiganol  prismatic  structure  with  Y-0  bond  ~  2.260  A  which  shortens  to  2.155  A  in  Y(acac)2. 
This  is  in  qualitative  agreement  with  the  structure  for  Y(dpm)3  obtained  from  gas  phase  electron 
diffraction  data  [32, 33].  The  effect  of  replacing  the  substituent  alkyl  group  on  the  ligand  (CH3  in 
acac  and  C(CH3)3  in  dpm)  with  H  atoms  (giving  the  3-oxo-propen-l-olate  ligand,  opo  ligand)  on 
the  ring  structure  was  investigated  and  found  to  be  negligible.  Hence  for  geometry  optimizations, 
the  terminal  group  on  the  ligand  was  replaced  with  H  atoms  to  determine  the  skeletal  ring 
stracture,  which  was  then  used  as  a  starting  guess  for  further  optimizations. 

The  acetylacetonate  radical  (acac)  can  exist  in  different  conformations  (3  in  the  keto  form 
and  1  in  the  enol  form)  [34].  The  different  conformations  are  shown  in  figure  2  and  the  relative 
energies  for  the  structures  are  given  in  table  1. 


Figure  2.  Different  possible  conformations  for  the  acac  ligand. 
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Conformation 

E° 

(hartrees) 

ZPE 

(hartrees) 

E°  +  ZPE 
(hartrees) 

Relative  Energy 
(kcal/mol) 

acac  1. 

-345.0892 

0.1082 

-344.9810 

0.0 

acac  2. 

-345.0737 

0.1092 

-344.9644 

10.4 

acac  3. 

-345.0933 

0.1091 

-344.9842 

-2.2 

acac  enol  4. 

-345.1065 

0.1099 

-344.9966 

-9.8 

Table  1.  Relative  energies  for  the  different  conformations  for  the  acac  ligand  radical. 


Precursor  Decomposition  Pathways 
1.  Copper  Precursors: 

Decomposition  of  Cu(acac)2  and  Cu(dpm)2  can  be  initiated  via  two  different  pathways: 

(a)  Dissociation  of  the  substituent  (CH3  for  acac  and  t-butyl  for  dpm)  on  the  ligand  ring  to  give 
the  corresponding  radical,  leaving  the  central  ring  intact. 

(b)  Dissociation  of  the  ligand  from  the  central  copper  atom  causing  the  breakup  of  the  ring 
structure  leaving  the  bare  metal  atom  behind.  This  may  occur  in  a  single  step  or  in  multiple 
steps. 

The  energy  required  to  break  the  necessary  bonds  gives  the  relative  strength  of  the  bonds 
and  the  bond  dissociation  energies  (BDE)  can  be  used  to  determine  the  initial  step  in  the  thermal 
decomposition  of  these  precursor  molecules. 


Figure  4.  Optimized  structure  for  Cu 
Figure  3.  Optimized  structure  for  (acac)2.  p-diketonate,  with  terminal  R  replaced  by  H. 


17 


The  nature  of  the  Cu-0  bond  does  not  change  with  the  change  in  the  terminal  alkyl  group 
on  the  ligand  and  the  bond  strength  of  the  Cu-O  bond  should  be  similar  for  the  different  Cu  P- 
diketonates.  In  order  to  reduce  the  computational  cost  involved  in  this  study,  the  effect  of 
replacing  the  terminal  t-butyl  groups  by  H  as  well  as  CH3  was  studied.  As  reported  earlier,  this 
has  no  influence  on  the  geometry  of  the  ring  structure  around  the  central  copper  atom. 

The  mean  Cu-0  strength  was  determined  as  follows 

BDE  (Cu-O)  =  [E(Cu.L2)  -  E(Cu)  -  2  E  (L)  ]  /  4 
where  E  is  the  zero-point  corrected  ground  state  energy. 


Figure  5.  Optimized  structure  for  Cu(dpm)2 


Species 

R 

BDE  (Cu-O) 
(NoZPEcorr.) 

BDE(Cu-O) 
(ZPE  corrected) 

BDE(Cu-O) 

Experimental 

[35] 

Cu(opo)2 

H 

5.2 

3.8 

- 

Cu(acac)2 

CH3 

40.1 

38.5 

39.9  ±2.4 

Cu(dpm)2 

C(CH3)3 

39.1 

37.9 

40.2  ±2.4 

Cu(hfac)2 

CF3 

42.8 

41.8 

39.2  ±2.6 

Table  2.  Mean  Bond  Dissociation  Energies  (BDE)  for  different  substituent  groups  on  the  Cu 
P-diketonate  precursor.  Energy  is  in  kcal/mol. 


The  mean  bond  dissociation  energy  for  Cu-0  bond  in  Cu(acac)2  and  Cu(dpm)2  is  ~  40 
kcal/mol.  The  calculated  BDE  energy  for  Cu(opo)2  gives  a  value  of  ~  4  kcal/mol.  Since  the 
terminal  alkyl  groups  do  not  influence  the  Cu-0  bond  strengths  considerably,  these  results 
indicate  that  the  substitution  of  the  t-butyl  group  by  an  H  atom  is  invalid.  The  reason  is  apparent 
when  one  considers  the  products  of  the  Cu-O  bonds  dissociation.  The  P-diketonate  bre^  up  to 
give  Cu  and  two  ligand  radicals.  The  radical  (CH3)3C-C(0)CHC(0)-C(CH3)3,  formed  from  the 


breakup  of  the  Cu(dpm)2,  is  a  stable  species,  whereas,  the  radical  HC(0)CHC(0)H,  is  highly 
unstable.  Calculations  done  with  R=CH3  give  excellent  agreement  with  the  energies  obtained 
with  R=C(CH3)3.  Hence  for  all  further  calculations,  the  t-butyl  can  be  replaced  by  CH3,  resulting 
in  considerable  savings  in  computational  time  and  effort. 

The  C-C  bond  strength  for  the  loss  of  the  t-butyl  group  from  Cu(dpm)2  is  72  kcal/mol. 
The  corresponding  bond  strength  for  the  removal  of  the  CH3  radical  is  90  kcal/mol.  These  values 
are  considerably  higher  than  the  mean  Cu-O  bond  strength  and  hence  Cu-0  bond  will  break  first 
followed  by  the  loss  of  the  t-butyl  or  the  methyl  group.  This  is  consistent  with  experimental 
observations  for  thermal  decomposition  of  Cu(dpm)2  where  the  IR  absorption  due  to  the  ring 
structure  begins  to  decrease  earlier  than  the  IR  absorption  for  the  C-H*  stretch  of  the  t-butyl  group 
[36, 37]. 

The  decomposition  is  thus  initiated  by  the  break-up  of  one  Cu-0  bond  to  form  a  complex  with 
one  bidentate  ligand  and  one  monodentate  ligand.  The  monodentate  ligand  rearranges  itself  to  a 
more  stable  conformation  and  the  resulting  complex  is  depicted  in 


Figure  6.  Structure  for  the  Cu(acac)2  with  one  monodentate  ligand  after  the  break-up  of  the  Cu-0 
bond. 


Fig.  7  (a) 


Fig.  7(b) 


Figure  7.  Optimized  structures  for  the  loss  of  CH3  radical  from  the  monodentate  ligand  of 
Cu(acac)2.  (a)  CH3  close  to  the  Cu  atom  is  lost  (b)  CH3  away  from  the  Cu  atom  is  lost. 
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The  second  step  is  the  splitting  of  the  CH3  radical  from  the  monodentate  ligand.  This  can 
take  place  in  two  different  ways:  (a)  CH3  close  to  the  central  Cu  atom  breaks  off  or  (b)  the  CH3 
further  from  the  Cu  atom  breaks  up.  Depending  on  which  CH3  splits  off,  we  may  get  different 
products. 

For  case  (a),  the  Cu-0  bond  between  the  monodentate  ligand  and  the  Cu  atom  breaks  and 
we  get  CH3(C0)CH(C0)  radical  and  Cu(acac)  which  can  further  break  up  to  give  acac  ligand  and 
Cu.  For  case  (b),  the  intact  bidentate  acac  ligand  can  abstract  an  H  from  the  monodentate  ligand 
forming  Hacac,  and  the  monodentate  ligand  breaks  up  forming  ketene  (CH2=C=0)  and  CH3CO 
radical.  The  initial  Cu-0  bond  splitting  is  the  rate-determining  step.  The  final  products  are  the 
Cu  atom  and  hydrocarbons.  The  chemistry  of  the  hydrocarbons  is  fairly  well  known  from 
combustion  chemistry  [12,  38].  The  bare  Cu  atom  gets  oxidized  to  give  CuO  and  CU2O. 

In  the  case  of  Cu(dpm)2,  we  get  t-butyl  radicals  (instead  of  CH3)  which  can  further 
disproportionate  to  give  butene  and  butane,  which  have  been  observed  experimentally  during 
therm^  decomposition  of  Cu(dpm)2  using  FTIR  spectroscopy  [39, 40]. 

The  above  observations  can  be  used  to  formulate  a  kinetic  mechanism. 


Copper  Acetvlacetonate: 


Cu(acac)2 

Cu(acac*)(acac) 

Cu(acac*)(acac) 

^  Cu(acac)  (CH3COCHCO)  +  CH 

Cu(acac)  (COCHCOCH3) 

^  Cu  +  CH2CO  +  Hacac  +  CHCO 

Cu(acac) -1- CH3COCHCO 

**  Cu(acac) 

Cu  +  acac 

Cu  +  O2 

->  CuO, 

Cu  +  CUO2 

->  2  CuO 

acac*  denotes  the  ligands  as  being  monodentate. 


CH2CO,  CHCO,  Hacac,  CH3COCHCO  are  hydrocarbons  whose  decomposition  and  oxidation 
chemistry  is  well  known  from  combustion  literature. 

Copper  Dipivaolvlmethanate: 


Cu(dpm)2 

-> 

Cu(dpm)2 

Cu(dpm)(COCHCOC4H9) 

C4H9COCHCO 

2C4H9 

**  Cu(dpm) 

2C4H9 

Cu  +  O2 

Cu  +  CUO2 

-> 

Cu(dpm)  -H  (COCHCOC4H9)  -(-  C4H9 
Cu(dpm)(COCHCOC4H9)  +  C4H9 
C  u  +  Hdpm  -I-  C4H8  +  CHCO  +  CO 
C4H8  -I-  CH2CO  +  CO 
C4H8  +  C4H10 
Cu  +  dpm 
C4H8  +  C4H10 
CUO2 
2  CuO 


**  Cu(acac)  and  Cu(dpm)  could  also  dissociate  via 


Cu(acac)  CuO  -1-  hydrocarbons 

Cu(dpm)  CuO  +  hydrocarbons 
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These  are  believed  to  be  minor  pathways  as  experiments  on  deposition  of  CuaO  from 
Cu(hfac)2  and  isotopically  labeled  H2**0  have  shown  that  the  oxygen  in  the  deposited  film  arises 
mainly  from  H2^*0  and  not  from  the  oxygen  in  the  ligand  molecule  [41]. 

Recently,  there  have  been  a  few  experimental  studies  to  determine  the  decomposition 
kinetics  for  the  metal  beta-diketonates.  Most  of  the  studies  have  been  based  on  mass 
spectrometry  where  the  concentration  of  the  observed  m/e  peaks  is  followed  to  fit  the  data  to  an 
Arrhenius  form  for  the  overall  decomposition  process.  The  kinetic  parameters  obtained  from 
these  studies  are  summarized  in  table  3. 


Precursor 

Environment 

Pre-exponential 

(1/sec) 

Activation  Energy 
(KcaPmol) 

Reference 

Cu(dpm)2 

vacuum 

1.5  X  10*' 

37.2  V  1.6 

[39] 

1.5  X  10^ 

27.0 

[42] 

9.6  X  10" 

32.7 

Cu(acac)2 

vacuum 

17.6 

[42] 

H2 

31.8 

02 

7.2  X  10' 

17.2 

Cu(acac)2 

vacuum 

I^EEQEHH 

27.5 

[43] 

Table  3.  Kinetic  parameters  of  the  gas-phase  pyrolysis  of  Cu  (P-diketonates). 


The  reported  values  are  observed  to  differ  widely  from  each  other.  Also,  the  kinetic 
parameters  are  based  on  mass-spectrometric  studies  may  not  be  able  to  capture  the  actual 
decomposition  process  which  is  believed  to  occur  via  a  radical  mechanism. 

In  light  of  the  conflicting  values  reported  in  the  literature,  quantum  chemistry  studies  can 
be  used  to  determine  the  kinetic  parameters.  The  rate-determining  step  for  both  Cu(acac)2  and 
Cu(dpm)2  decomposition  appears  to  be  the  Cu-0  bond  break-up.  This  step  is  assumed  to  occur 
without  a  transition  state.  Based  on  the  quantum  chemistry  calculations  for  the  copper 
diketonates  with  the  bidentate  and  the  monodentate  ligands,  the  activation  energy  for  this  step  is  ~ 
27.0  kcal/mol.  Pre-exponentials  for  this  type  of  reactions  range  from  10^^  to  10*“*.  The  value  of 
the  pre-exponential  of  10^^  seems  reasonable.  These  parameters  can  be  used  in  the  finite  element 
simulations  to  model  the  CVD  deposition  of  the  HTS  thin  films. 

Decomposition  in  the  Presence  of  Oxygen: 

The  effect  of  oxidant  gas  such  as  oxygen  on  the  decomposition  of  the  Cu  p-diketonates 
was  also  investigated.  Oxygen  binds  to  the  Cu(opo)2  precursor,  forming  a  hexa-coordinated 
complex  as  shown  in  Figure  4.  This  is  accompanied  by  the  weakening  of  the  Cu-0  bond  between 
the  central  copper  atom  and  the  oxygen  atom  from  the  ligand.  The  bond  length  between  Cu-0 
increases  from  1.94  A  in  Cu(opo)2  to  2.14  A  in  Cu(opo)2.02,  which  is  indicative  of  the  weaker 
Cu-O  (ligand)  bond.  Cu(opo)2  has  the  H  atoms  as  substituents  on  the  ligand.  It  is  possible  that 
such  a  structure  is  unstable  and  falls  apart  easily.  It  seems  that  oxygen  reacts  with  the  metal 
directly  forming  Cu02  and  CuO. 

In  several  deposition  studies,  where  Cu(hfac)2  is  employed  as  a  copper  precursor,  water 
vapor  is  used  along  with  O2  as  the  oxidant.  Quantum  chemistry  calculations  show  the  tendency 
of  a  Cu(hfac)2  to  bind  to  the  oxygen  of  the  water  to  form  a  penta-coordinated  complex.  This 
tendency  of  flourinated  copper  precursors  to  coordinate  an  nucleophillic  nucleus,  such  as  O  in 
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H2O,  to  form  higher  coordinated  compounds  has  also  been  observed  experimentally,  where  a 
water  molecule  was  shown  to  bind  to  the  Cu  atom  of  Cu(hfac)2  [41,  44].  Recent  experiments  on 
deposition  of  CU2O  from  Cu(hfac)2  and  isotopically  labeled  H2**0  have  shown  that  the  oxygen  in 
the  deposited  film  arises  mainly  from  H2**0  and  not  from  the  oxygen  in  the  ligand  molecule  [41]. 


Figure  8.  Optimized  structure  of  the  Cu(opo)2-02.  Coordination  with  O2  increases  the  Cu-O 
(ligand)  bond  distance. 


Figure  9.  Optimized  structure  of  the  Cu(hfac)2-H20. 


2.  Yttrium  Precursors: 

Yttrium  tris(dipivaloylmethante)  or  Y(dpm)3  is  the  most  common  precursor  employed  for 
deposition  of  yttrium  containing  HTS  films.  The  large  number  of  electrons  present  in  Y  along 
with  the  bulky  r-butyl  ligands  makes  the  task  of  studying  this  molecule  quite  challenging. 
ONIOM  calculations  were  done  to  obtain  the  structure  of  Y(dpm)3  and  compared  with  the 
structure  for  Y(acac)3,  where  the  t-butyl  groups  have  been  replaced  by  CH3. 
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Figure  10.  Optimized  structure  of  the  Y(acac)3. 

Y(acac)3  was  chosen  as  a  model  for  Y(dpm)3  and  bond  dissociation  energies  were 
calculated  for  the  Y-0  bond.  The  bond  dissociation  energy  for  the  removal  of  the  first  ligand 
from  Y(acac)3  according  to  the  reaction 

Y(acac)3 Y(acac)2  +  acac 

is  61  kcal/mol.  The  mean  bond  dissociation  energy  for  the  Y-0  in  Y(dpm)3  is  68  kcal/mol.  The 
initial  step  is  the  decomposition  of  Y(dpm)3  is  the  splitting  of  the  bulky  t-butyl  ligand.  The  loss 
of  t-butyl  from  the  ligand  causes  the  Y-0  bond  to  weaken.  This  weakening  of  the  bond  causes 
the  ligand  to  be  split-off  from  the  molecule  resulting  in  Y(dpm)2  which  is  more  stable  compared 
to  the  parent  molecule,  Y(dpm)3.  This  is  evidenced  from  the  shortening  of  the  Y-O  bond  in 
Y(dpm)2  and  Y(acac)2  as  compared  to  the  Y(dpm)3  and  Y(acac)3.  Y(dpm)2  has  been  detected  as  a 
stable  species  in  the  gas  phase  at  temperatures  as  high  as  760  K  [33, 45] . 

Two  different  structures  were  obtained  for  Y(acac)2.  In  structure  I,  Y  is  square-planar 
coordinated,  and  in  structure  H,  Y  is  in  tetrahedral  coordination.  The  two  structures  are  shown  in 
figure  11.  A  planar  structure  has  also  been  predicted  for  Y(dpm)2  based  on  electron  diffraction 
data  [33].  Quantum  chemistry  calculations  predict  the  tetrahedral  stmcture  to  be  more  stable 
compared  to  the  planar  structure  by  38  kcal/mol. 


Figure  11.  Optimized  structures  for  Y(acac)2. 


A  kinetic  mechanism  for  decomposition  of  Y(dpm)3  based  on  the  above  observations  is: 

Y(dpm)3 

->  (Y(dpm)3-C4H9)  +  C4H9 

Y(dpm)3 

^  ( Y(dpm)3  -  3  C4H9  )  +  3  C4H9 

(Y(dpm)3-C4H9) 

^  Y(dpm)2  +  C4H9+(C0)CH(C0) 

Y(dpm)2 

Y  +  2dpm 

4Y  +  302 

2Y2O3 

Conclusion  on  CVD  Kinetics  Modeling 

We  have  developed  a  systematic  methodology  to  look  at  the  structure  and  reactivity  of  the  metal 
P-diketonates  precursors,  used  in  the  chemical  vapor  deposition  of  high  temperature 
superconductors,  using  first  principle  quantum  chemistry  methods.  The  large  size  of  the 
precursor  molecules  has  so  far  limited  studies  of  this  nature.  Very  little  information  is  available 
on  the  chemistry  of  these  compounds.  Semi-empirical  methods  were  found  inadequate  in 
describing  the  structures  and  traditional  ab-initio  methods  were  deemed  to  be  very  expensive, 
given  the  size  of  the  molecules.  We  have  tested  and  validated  different  density  functionals  for 
describing  these  molecules.  Based  on  our  validations,  an  effective  core  potential  was  chosen  to 
describe  the  metal  (Cu  and  Y).  We  have  investigated  the  decomposition  mechanism  for 
Cu(acac)2  and  Cu(dpm)2,  two  of  the  most  common  precursors  for  copper  in  HTS  films.  A  kinetic 
mechanism  for  decomposition  was  formulated  and  the  initial  break-up  of  the  Cu-O  bond  was 
determined  to  be  the  rate-limiting  step.  A  similar  investigation  was  also  carried  out  for  Y (acac)2, 
which  serves  as  a  model  for  Y(dpm)3.  A  preliminary  mechanism  has  been  put  together  based  on 
the  work  done  so  far. 
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3.  Reactor  Scale  Modeling  of  MOCVD  of  HTS 

The  objective  of  this  work  is  to  develop  a  detailed  and  accurate  reactor-scale,  physical  model  that 
predicts  deposition  rate  and  stoichiometry  along  the  wafer  surface  when  the  operating  conditions 
and  inflow  gas  composition  are  specified.  The  physical  model  has  to  address  precursor 
decomposition  and  oxide  formation  chemistry,  and  the  transport  processes  (fluid  mechanics, 
species  diffusion,  and  heat  transfer)  associated  with  CVD.  The  reactor  modeled  here  is  the 
Thomas  Swan  reactor  used  by  STI  for  their  production  of  YBCO  superconducting  thin  films 
(HTS).  For  modeling  chemistry,  we  have  used  kinetics  data  developed  from  quantum  chemistry 
calculations  described  earlier,  supplemented  by  other  data  in  the  literature. 


A  detailed  description  of  various  aspects  of  non-metal  CVD  is  the  subject  of  a  recent  book  [59]. 
Various  reactor  configurations  such  as  planetary,  barrel,  and  pancake  reactors  for  MOCVD  are 
described  in  [60].  The  more  specific  topic  of  MOCVD  of  high-Tc  HTS  is  discussed  in  two 
excellent  review  papers  [61,  62],  and  there  are  several  reports  of  experimental  results  in  the 
literature,  e.g.,  [67,68]. 

MOCVD  of  HTS  materials  is  a  relatively  new  area  of  research,  and  there  is  little  published  on 
reactor-scale  simulation  for  comparison  with  our  work.  Hence,  a  step-by-step  approach  towards 
acquiring  capability  for  reliable  MOCVD  simulations  was  adopted  in  this  study.  The  first  step 
was  to  perform  simple  CVD  calculations  with  one-step,  finite-rate  chemistry  for  the  above  2-D 
reactor  and  compare  with  results  in  the  literature.  The  case  studied  was  silicon  epitaxy  using 
trichlorosilane  (SiHCls),  and  the  detailed  problem  was  identical  to  that  solved  by  the  Habuka  et  al 
[52].  A  commercial  software  package,  CFD-ACE  [63],  popular  in  the  semiconductor  industry, 
was  used  for  our  modeling  work.  The  equations  solved  by  CFD-ACE  to  calculate  CVD 
deposition  rates  are  described  next  using  the  silicon  epitaxy  example,  followed  by  the  detailed 
validation  problem  before  describing  MOCVD  of  HTS. 


Equations  of  CVD 

The  global  CVD  reaction  for  deposition  of  epitaxial  silicon  at  the  surface  of  the  wafer  is 
represented  by: 

SiHCls  +  H2  Si  +  3HC1 
(gas)  (gas)  (solid)  (gas) 


An  Arhennius-type  kinetic  rate  constant,  k,  for  the  global  one-step  reaction  was  obtained  from 
Habuka  et  al  [5]  as: 


k  =  2650exp 


-138,000 

RT 


;  (units :  m*  l(gm  -  mole.s)) 


Here,  Tis  the  temperature  in  Kelvin,  ka  is  the  Boltzman  constant.  The  deposition  rate  is  then 
obtained  using 


J=Jk[SiHCl3]  [h] 

Here,  the  quantities  in  square  brackets  are  molar  concentrations,  d  is  the  deposition  rate  in 
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moles/m  Vs.  The  equations  that  are  solved  using  CFD-ACE  are  the  mass,  momentum,  energy,  and 
species  concentration  equations.  The  mass  conservation  (continuity)  equation  has  the  following 
form: 


(/TUj)  =  0 


The  mixture  density,  p,  depends  on  temperature,  and  is  computed  using  the  ideal  gas  law.  The 
first  term  in  this  and  the  following  equations  represent  time-dependence  of  the  quantities. 
Although  all  the  flow  and  temperature  fields  are  steady-state,  the  term  is  included  here  for 
completeness.  The  term  Uj  is  the  f'  component  of  the  instantaneous  velocity.  The  momentum 
conservation  equations,  or  the  Navier-Stokes  equations  have  the  following  form: 


dt 


(pu.) 


+ 


dp  ^  d  dui  ^  duj  2  ^ 

3Xi  3Xj.  ^  dxj  dx,  3^  dx^. 


+  f>f, 


In  the  above  equation,  p  is  the  static  pressure,/-  is  the  body  force,  and  is  the  Kronecker-delta 
function.  The  mixture  viscosity,  p,  is  also  temperature  dependent.  The  three  momentum  equations 
yield  the  three  Cartesian  components  of  velocity.  The  energy  balance  equation  is  written  in  terms 
of  the  enthalpy,  h  (static  enthalpy  for  these  relatively  low  speed  flows),  as  follows: 


dx- 


dXj 


The  last  term  on  the  right  side  represents  viscous  dissipation,  and  is  not  considered  in  our 
calculations  because  it  is  very  small  at  these  low  velocities.  The  second  term  on  the  left  hand  side 
is  the  enthalpy  convection  term.  The  first  term  on  the  right  hand  side  represents  conduction  heat 
flux,  which  is  modeled  using  Fourier’s  law.  Replacing  the  temperature  variable,  T,  by  h  using  h 
=  Cp  T,  where  Cp  is  the  specific  heat  capacity,  this  term  takes  the  following  form: 

kr  dT 

^  — _ _ 


Here,  kj-is  the  thermal  conductivity  of  the  gas.  The  species  transport  equation  has  the  following 
form 


The  term  is  the  mass  fraction  of  the  i*  species,  J,-  is  its  diffusive  mass  flux,  and  Si  is  a 
volumetric  source  term  representing  creation  or  consumption  of  the  i*  species  via  chemical 
reaction.  Since  the  CVD  problem  here  is  solved  as  a  surface  reaction,  S,  =0.  The  above  equation 
is  solved  for  any  two  of  the  three  gaseous  species,  and  the  third  mass  fraction  is  obtained  by 
subtracting  the  sum  of  the  rest  from  unity.  The  diffusive  mass  flux,  /„  is  the  sum  of  the  flux 
driven  by  concentration  gradient ,  7,-*^,  and  the  thermodiffusive  (Soret)  flux,  7,-^  driven  by 
temperature  gradient.  The  concentration  diffusive  flux  for  any  species  depends  on  the  mass 
concentration  gradients  of  all  the  species  through  the  Stefan-Maxwell  equations  [67]  shown 
below: 
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dcOj 

dxj 


+  0)i- — (inAf)  = 


N 

y 

k*i 


M  (OJI 


M. 


M; 


Here,  Dik,  is  the  binary  diffusion  coefficient  of  species  i  with  respect  to  species  k.  M-,  is  the 

molecular  weight  of  species  i,  M  is  the  average  molecular  weight  of  the  gas,  and  N  is  the  total 
number  of  species.  The  above  equation  is  solved  iteratively  to  yield  jf  as  a  function  of  the 
various  mass  fractions  subject  to  the  following  constraint. 

N 

=  0 
i=i 

The  Soret  diffusion  flux,  which  drives  the  gas  away  from  hot  walls  towards  cooler  walls,  is  given 
by 

Jj  D!  A(inr) 

Here,  is  the  thermodiffusive  coefficient.  As  stated  earlier,  the  temporal  terms  in  all  the  above 
equations  are  neglected  in  the  calculations.  This  pseudo-steady-state  condition  is  justified 
because  the  time  scales  of  transport  are  much  shorter  (of  the  order  of  a  few  seconds)  than  the  time 
scales  for  film  growth,  especially  at  lower  pressures.  The  boundary  condition  for  the  species 
equation  at  the  deposition  surface  is  obtained  by  balancing  the  sum  of  the  advective  and  diffusive 
fluxes  of  the  species  normal  to  the  wall  against  the  flux  of  species  generated  (or  consumed): 


P(po>,)+J,)n  +  =  0 

l=l 


Here,  V  is  the  velocity  vector,  n  is  the  normal  to  the  surface,  nr  is  the  total  number  of  surface 
reactions,  and  dn  is  the  molar  flux  of  species  i  generated  (or  consumed)  in  reaction  1. 

The  above  transport  equations  are  solved  using  a  control- volume  approach  [63]  for  the  structured 
grid.  In  this  solution,  a  first-order  upwind  scheme  was  used  for  the  convective  fluxes.  CFD-ACE 
uses  a  pressure-based  approach  to  solving  the  Navier-Stokes  equations  in  which  the  continuity 
equation  is  used  to  recast  these  three  equations  is  terms  of  pressure.  This  solution  used  the 
SIMPLEC  method,  a  variation  of  the  well-established  SIMPLE  algorithm  [64]. 


CVD  validation  studies  for  silicon  epitaxy 

Figure  12  shows  the  geometry  used  in  the  2-D  model  obtained  from  Habuka  et  al  [52].  The  end- 
to-end  length  of  the  radiantly-heated  chamber  is  0.705  m  and  the  total  height  is  0.4  m.  The  top 
part  of  the  left  wall  and  the  entire  right  wall  are  at  300  K,  and  the  sections  of  the  top  and  the 
bottom  wall  immediately  above  the  susceptor  are  at  specified  temperatures.  The  rest  of  the  walls 
are  adiabatic.  The  flow  enters  at  atmospheric  pressure  from  the  top  left  comer  of  the  chamber, 
and  exits  from  the  bottom  right  comer.  The  8”  wafer  is  located  at  the  center  of  the  12”  susceptor, 
0.205  m  from  the  left  end  of  the  chamber. 
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Figure  12.  Schematic  of  a  two-dimensional  horizontal  CVD  reactor  [52].  The  figure  is  not  to 
scale.  Shaded  walls  are  adiabatic.  The  part  of  the  walls  next  to  the  IR  heaters  are  semi-transparent 
and  at  elevated  temperatures. 


A  gaseous  mixture,  consisting  of  trichlorosilane  (SiHCls,  nominal  mass  fraction=0.71)  and 
hydrogen,  is  injected  into  the  chamber  at  room  temperamre  (300K)  and  a  velocity  of  0.67  m/s. 
The  wafer  temperature  is  isothermally  elevated  to  a  nominal  value  of  1423 K.  The  temperatures 
of  the  hot  sections  of  the  top  and  bottom  walls,  Twaii,  were  measured  by  Habuka  et  al  [52],  and 
expressed  as  the  following  linear  function  of  the  susceptor  temperature,  Tjus: 

Twaii  =  730  +  (770-730)(Tsus-1393)/(1453-1393) 

For  Tsus  =  1423K,  Twaii  =  750K.  The  other  sections  on  the  top  and  the  bottom  walls  are  kept  at 
300K.  The  susceptor  (and  hence,  the  wafer)  temperatures  are  assumed  to  be  independendy 
controlled  to  excellent  uniformity. 
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Figure  13.  Mesh  for  numerical  solution  using  CFD-ACE.  The  90  X  53  mesh  (with  3399  cells  or 
control  volumes  in  all)  is  clustered  near  the  walls  and  near  regions  of  high  temperature  gradients 
at  the  edges  of  the  wafer.  For  clarity,  the  vertical  dimensions  are  magnified  five  times  the 
horizontal  dimensions. 

Figure  13  shows  the  mesh  generated  for  the  control-volume  solution.  The  solution  converged 
after  600  iterations  in  about  half-an-hour  on  a  Pentium  PC.  Figure  14  shows  the  velocity  vectors 
with  superposed  temperatures  for  nominal  conditions  (trichlor  mass  fraction  of  0.71,  wafer 
temperature  of  1423K,  and  top  and  bottom  quartz  wall  temperatures  of  750K).  The  gas  is  heated 
up  considerably  by  the  susceptor  and  the  wall,  and  speeds  up  along  the  wafer  surface.  Figure  15 
shows  comparison  of  deposition  rate  uniformity  with  Habuka  et  al  [52].  Figure  16  shows  that 
wafer  rotation  significantly  improves  deposition  uniformity,  assuming  that  the  rotation  period  is 
much  smaller  than  the  deposition  period  (which  is  almost  always  the  case).  The  CFD-ACE 
results  compares  quite  well  with  those  Habuka  et  al  [52],  the  average  deposition  rate  difference 
being  about  10%. 
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Figure  15.  Deposition  profile  along  flow  direction.  Comparison  with  Habuka,  et  al  [52] 


Figure  16.  Effect  of  wafer  rotation  on  deposition  rate  uniformity. 

The  deposition  rates  were  calculated  by  varying  the  trichlor  mass  fraction.  The  results  are  plotted 
in  Figure  17.  These  results  agree  well  with  Habuka’s  experimental  data.  The  slight  over¬ 
prediction  of  the  deposition  rate  may  be  attributed  be  due  to  two  possible  causes.  First,  the 
temperatures  at  the  adiabatic  parts  of  the  wall  are  relatively  high,  and  in  reality,  there  is  probably 
some  deposition  on  the  walls  leading  to  reactant  depletion  downstream.  Second,  there  is  always  a 
small  amount  of  HCL  etching  of  deposited  silicon  that  reduces  the  overall  deposition  rate. 


Figure  17.  Average  deposition  rate  as  function  of  average  molecular  weight  of  gas  at  inlet. 
Nominal  molecular  weight  is  6.67  kg/kg-mole  (when  trichlor  mass  fraction  =  0.71). 


Several  tests  were  carried  out  on  the  model  to  establish  convergence  on  the  basis  of  mesh 
refinement  and  number  of  iterations  needed.  The  results  are  shown  in  Figure  18.  We  find  that 
reducing  iterations  from  5000  to  KXX)  results  in  an  average  difference  of  0.5%  in  the  deposition 
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rate.  Based  on  this  result,  it  is  arbitrarily  decided  that  for  this  geometry,  1500  iterations  are 
sufficient.  Reducing  the  number  of  cells  from  3399  to  1381  changed  the  average  deposition  rate 
by  0.8%,  and  the  average  horizontal  centerline  temperature  by  0.7%.  Hence,  1381  cells  are 
deemed  sufficient.  It  is  also  noted  that  the  time  per  iteration  rises  non-linearly  for  more  refined 
mesh,  which  is  to  be  expected  from  this  solver. 

From  this  validation  study,  we  concluded  that  the  results  for  a  simple  CVD  model  agree  well  with 
those  published  in  the  literature,  giving  us  the  necessary  confidence  to  proceed  with  modeling  of 
MOCVD  of  YBCO  thin  films. 


Figure  18  Convergence  study  using  deposition  rates  and  gas  temperatures. 
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Physical  Model  for  MOCVD  of  HTS 


A  2-D  version  of  the  geometry  of  the  commercial,  hot  wall  reactor  (Thomas  Swan)  used  at  STI 
was  developed  for  modeling  MOCVD  of  HTS.  The  reactor  is  infinite  in  the  transverse  direction 
and  is  shown  in  Figure  19.  The  center  of  the  wafer  of  diameter  2”  is  located  20  cm  from  the  left 
end  of  the  reactor  through  which  the  gas  mixture  enters.  The  ceiling  of  the  reactor  slopes 
downward  in  order  to  increase  species  concentration  gradients  in  the  vertical  direction,  which  in 
turn  increases  species  flux  towards  the  wafer  surface.  This  design  thus  improves  deposition 
uniformity  by  compensating  for  species  depletion  along  the  wafer  caused  by  deposition. 


4  cm 
10  cm 


20  cm 
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30  cm 


Figure  19.  Schematic  of  STFs  Thomas  Swan  MOCVD  reactor. 

Figure  20  shows  the  structured  mesh  constructed  for  the  above  reactor,  magnified  five  times  in 
the  vertical  direction  for. clarity.  The  mesh  consists  of  574  (14x41)  cells  which  was  found 
sufficient  for  obtaining  a  converged  solution. 


Figure  20.  Mesh  for  numerical  solution  using  CFD-ACE,  clustered  near  the  walls.  For  clarity, 
the  vertical  dimensions  are  magnified  five  times  the  horizontal  dimensions. 


For  chemistry,  we  have  used  one-step  global  mechanisms  for  precursor  decomposition.  For 
Cu(dpm)3,  the  kinetic  rate  constants  were  obtained  from  quantum-mechanical  calculations 
performed  in  this  program.  For  the  other  two  precursors,  the  kinetic  mechanism  were  obtained 
from  the  Colorado  School  of  Mines  [53],  another  VIP  program  participant.  They  derived  the 
'global’  empirical  mechanism  from  existing  literature.  The  overall  mechanism  consists  of  seven 
gas  phase  reactions; 


Y(dpm)3  +  O2  — ^ 


Y  +  3(dpm)  +O2 


Y(dpm)3  Y  +  3  (dpm) 
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Cu(dpm)2  — >  Cu  +  2  (dpm) 

Ba(dpm)2  — >  Ba  +  2  (dpm) 

4Y  +  3O2  ^  2Y2O3 
2  Ba  +  O2  — ^  2BaO 
2Cu  +  02-^  2CuO 

The  last  three  oxidation  reactions  are  relatively  fast.  The  Arrhenius  rate  constants  for  the  first 
four  reactions  are  as  follows.  The  pre-exponential  terms  for  the  four  reactions  are  5.9xlO  '^ 
m^/kg-mole.s,  2.2x10^  s^  l.OxlO'^  s'\  and  2.2x10^  s'\  respectively.  The  characteristic 
temperature  (EJR)  are  13590.4  K,  17264.9  K,  13620.6  K,  and  17264.9  K,  respectively. 

The  oxides  produced  in  the  gas  phase  reaction  diffuse  to  the  surface  and  their  deposition  was 
modeled  using  sticking  coefficients.  Since  the  transport  is  diffusion  limited,  choice  of  sticking 
coefficient  is  not  critical.  This  issue  is  discussed  later  in  the  section.  Sticking  coefficients  of 
unity  were  used  for  all  three  oxide  species.  The  transport  properties  of  some  of  the  species  were 
calculated  using  Lennard-Jones  parameters.  For  the  oxides,  data  on  melting  point  temperatures, 
Tm  [65],  and  the  liquid  molar  volume,  Vb,  were  used  to  calculate  the  parameters  using  the 
following  correlations  from  the  literature  [65]. 

e/k  =  1.92T,„;  0=1.166  Vb‘" 

For  the  precursors,  melting  point  data  was  used  to  calculate  e/k.  The  collision  diameter,  o,  was 
calculated  from  data  on  Oy  for  precursor  and  oxygen  [66].  Since  data  for  Cu(acac)2  is  reported, 
we  use  these  values  for  Cu(dpm)2  for  lack  of  more  accurate  values  of  the  transport  properties. 
Finally,  the  values  of  Lennard-Jones  parameters  used  in  the  model  are  tabulated  below. 


Species 

Melting 

pt(K) 

S/ks 

(K) 

a 

(A) 

Y(thd)3 

523 

601 

17.7 

Ba(thd)2 

606 

697 

15.3 

Cu(thd)2 

625 

719 

10.8 

thd 

111 

9.38 

Y203 

2410 

46in 

4.35 

BaO 

1918 

3683 

3.66 

CuO 

1326 

2546 

2.83 

The  oxidation  reactions  are  instantaneous  while  the  decomposition  reactions  are  finite-rate,  make 
the  system  of  equations  very  stiff.  We  decided  to  simplify  the  model  by  eliminating  the  four 
oxidation  reactions.  The  metallic  atoms  are  assigned  the  transport  properties  of  the  oxides,  which 
then  diffuse  to  the  surface.  Because  the  mole  fractions  of  the  precursors  are  very  small  (of  the 
order  of  tens  of  ppm),  the  excess  oxygen  not  used  up  in  oxide  formation  are  also  very  small,  and 
its  effect  on  the  flow  and  temperature  fields  is  insignificant.  This  simplification  halved  the 
computational  time  without  any  loss  in  accuracy.  A  further  saving  of  25%  in  simulation  time  was 
achieved  by  neglecting  Soret  diffusion  (thermodiffusion)  and  Stefan-Maxwell  conservation. 


34 


Because  the  temperature  gradients  inside  the  reactor  are  small,  Soret  diffusion  was  found  to  be 
negligible.  The  concentrations  of  the  three  oxides  are  very  small.  Hence,  neglecting  Stefan- 
Maxwell  conservation  for  multi-species  diffusion  caused  negligible  error  in  the  deposition  rate. 
A  typical  simulation  would  converge  within  1200  iterations  in  approximately  two  hours  of  CPU 
time  on  Pentium  PC. 

Results  are  shown  below  for  the  following  representative  steady-state  operating  conditions  [53]. 
The  precursors  are  well  mixed  with  oxygen,  nitrogen,  and  argon  and  enter  the  reactor  at  10  Torr 
at  a  velocity  of  2  m/s  and  temperature  of  513K.  The  inlet  mass  fractions  are:  02=0.456, 
N2=0.422,  Ar=0.114,  Y(dpm)3=9.32xl0"‘,  Ba(dpm)2=1.57xl0■^  Cu(dpm)2=5.02xl0'l  The  wall 
of  the  plenum  is  insulated  up  to  the  highest  point  of  the  chamber  (10  cm  from  the  entry  port).  The 
remaining  20  cm  of  the  reactor  are  heated  radiantly  so  that  both  the  top  and  bottom  surfaces  are 
kept  at  1073K.  Although  the  2”  diameter  wafer  is  centered  20  cm  from  the  entry  port,  deposition 
occurs  along  the  entire  length  of  the  heated  lower  surface.  Figure  21  shows  the  flow  velocity 
vectors  with  the  color  scheme  showing  the  magnitude  (speed).  The  flow  speeds  up  considerably 
as  it  gets  hotter.  Figure  22  shows  the  temperature  contours,  and  Figure  23  shows  the  Y2O3  mass 
fraction  contours. 

In  this  model,  YBCO  deposition  rate  is  calculated  by  adding  the  deposition  rates  of  each  of  the 
three  individual  oxides. 
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Figure  23.  Yttrium  oxide  mass  fraction  distribution  in  reactor. 


Figure  24  shows  the  precursor  concentrations  along  the  reactor  length  5  mm  above  the  wafer 
surface,  and  also  the  oxide  concentrations  across  the  height  reactor  at  the  center  of  the  wafer.  It  is 
seen  that  the  precursors  decompose  very  quickly  as  they  enter  the  heated  region. 


Figure  24.  Spatial  profiles  of  various  species  along  reactor  length. 
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Figure  25.  Spatial  profiles  of  various  species  in  the  vertical  direction. 

Figure  25  shows  the  oxide  mass  fractions  in  the  vertical  at  the  center  of  the  wafer.  The  oxides 
concentrations  are  very  small  near  the  surface  because  of  depletion  due  to  deposition. 


Figure  26.  YBCO  deposition  rates  in  Angstroms/minute,  wafer  not  being  rotated. 
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Figure  27.  YBCO  deposition  rates  in  Angstroms/minute,  wafer  rotated. 


Figure  26  shows  the  deposition  rate  along  the  diameter  of  the  2”  wafer  for  the  case  where  the 
wier  is  static,  while  Figure  27  shows  deposition  uniformity  when  the  wafer  is  rotated.  Although 
the  deposition  occurs  more  uniformly  with  wafer  rotation,  the  20%  non-uniformity  within  a  2” 
wafer  is  still  rather  high.  The  average  deposition  rates  of  35.7  A/min  are  in  the  range  of  30-40 
A/min  deposition  rate  obtained  by  STL 

Earlier  it  was  stated  that  the  deposition  process  is  transport  limited.  Figure  28  shows  that  this  is 
indeed  the  case,  with  deposition  rate  of  yttrium  oxide  showing  little  change  with  a  five-fold 
increase  in  sticking  coefficient.  Deposition  on  the  entire  heated  substrate,  rather  than  just  the  2” 
wafer,  is  shown  in  this  graph.  Figure  29  shows  the  stoichiometry  of  YBCO  thin  film  across  the 
wafer.  The  stoichiometry  changes  fi'om  point-to-point.  In  general,  its  is  desirable  to  have  the 
BaO/YjOs  ratio  slightly  below  two,  and  the  CuOA^’aOs  ratio  slighdy  above  three  [62]. 
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4.  POD-Galerkin  for  Chemistry  Reduction 


The  motivation  for  this  research  stems  from  the  need  to  develop  reduced  models  of  chemistry  in 
reacting  flow  problems.  The  numerical  solution  of  such  problems,  where  the  right-hand-side  is 
complicated  by  the  presence  of  homogeneous  chemical  reaction,  can  be  difficult  and 
computationally  expensive.  Convergence  of  the  numerical  solutions  can  become  even  stiffer  in 
the  case  of  heterogeneous  chemical  reactions,  which  introduce  highly  nonlinear  boundary 
conditions  to  the  sets  of  nonlinear  PDEs  whose  solution  is  desired;  examples  include  catalytic 
combustion  as  well  as  crystal  growth  by  CVD.  The  complexity  of  the  chemistry  clearly  increases 
with  the  number  of  gas-phase  and  surface  species  involved  in  the  reactions.  It  is  therefore 
understandable  that  in  reacting  flows  the  reduction  of  the  chemistry  to  a  few  components 
(pseudospecies)  will  not  only  improve  the  speed  of  the  solution  by  decreasing  the  number  of 
equations  to  be  solved,  but  will  also,  especially  in  the  case  of  complex  heterogeneous  reactions, 
dramatically  improve  convergence. 

The  Proper  Orthogonal  Decomposition  (POD)  method  was  used  for  chemistry  reduction.  While 
PODs  have  been  used  before  for  reducing  the  order  of  various  transport  fields  (temperature, 
velocity),  the  parallel  individual  reduction  of  chemistry  was  the  issue  we  needed  to  investigate.  A 
comparison  with  the  simultaneous  construction  of  PODs  from  all  transport  and  chemistry  fields 
needs  to  be  also  constructed  in  also  to  assess  the  effectiveness  of  the  technique. 


POD-based  compression  of  two  kinetic  schemes  were  considered;  the  first,  describing  the 
MOCVD  growth  of  Si  from  SiClaHa  using  H2  as  the  carrier  gas,  and  the  second  describing  Si 
deposition  kinetic  mechanism  in  the  SiC^Ha  /  B2H6  /  H2  system.  As  a  first  step  we  examined  the 
POD  content  of  the  concentration  dynamics  and  of  the  temperature  dependence  of  a  subset  of  the 
first  reaction  mechanism.  For  this  purpose  the  gas-phase  reaction  mechanism  was  considered 
separately.  This  scheme  involves  5  gaseous  species  and  6  chemical  reactions;  namely: 

1.  SiCl2H2  ~>  SiCl2  +  H2 

2.  SiCl2  +  H2  ~>  SiCl2H2 

3.  SiCl2+H2~>  HSiCl  +  HCl 

4.  HSiCl  +  HCl  ->  SiCl2  +  H2 

5.SiCl2H2~>HSiCl  +  HCl 

6.  HSiCl  +  HCl  ->  SiCl2H2 

An  “effective”  batch  reactor  chemistry  was  taken  into  account  and  the  corresponding  system  of 
ode's  was  constructed.  The  concentration  curves  of  SiCl2H2  from  the  time-integration  are 
shown  in  Figure  30  (solid  lines)  for  3  different  temperatures  T=950,1100  and  1200  K  and  an 
operating  pressure  of  1  atm.  A  data  ensemble  of  time-dependent  concentration  curves  of  the  5 
gaseous  species  for  the  3  different  temperatures  was  considered  and  the  eigenvalues  of  the  two- 
point  correlation  matrix  revealed  that  one  mode  captures  99.99%  of  the  system's  energy.  The 
projections  of  the  transient  evolution  of  SiCl2H2  concentration  (obtained  from  the  ODE  model)  on 
the  corresponding  eigenfunction  are  depicted  in  Figure  30  (dashed  lines).  As  can  be  seen  both 
the  transient  concentrations  and  the  steady  states  at  different  temperatures  are  accurately 
captured.  This  is  an  indication  that  the  kinetic  mechanism  of  the  Si  deposition  system  can  be 
successfully  reduced  via  the  POD  method. 
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Initial  results  for  the  Galerkin  projection  of  these  ODEs  on  the  POD  modes  indicate  that  the 
short-term  dynamics  are  accurately  captured,  but  show  deviations  for  the  long-term  dynamics. 
As  the  research  concludes,  we  are  studying  the  possibility  of  creating  a  richer  data  ensemble, 
perturbing  the  data  in  directions  orthogonal  to  the  most  important  POD  vectors. 


Figure  30.  Concentration  transients  of  SiCl2H2  obtained  from  the  ODE  system  (solid  lines)  and 
reconstruction  of  the  concentration  transients  from  the  projections  of  the  data  on  the  POD  modes 
(dashed  lines  for  three  different  temperatures  T=950, 1100  and  1200  K. 
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5.  Model-Based  Control  for  HTS  MOCVD 

This  section  contains  ideas  for  real-time  and  run-to-run  controller  design  of  a  CVD  reactor. 
Physical  Model 

A  2D  model  has  been  derived  that  approximates  the  static  global  behavior  of  the  commercially 
available  STI  reactor,  see  references  [25,26].  This  model  constitutes  seven  “global”  gas  phase 
reactions,  a  mass  conservation  equation,  momentum  conservation  equations,  an  energy  balance 
equation,  concentration  diffusive  flux  equations  (these  are  all  pdes)  and  several  algebraic 
relations.  Most  of  these  equations  are  strongly  coupled.  So  far,  the  dynamical  part  of  these 
equations  has  been  neglected. 

Control  Objectives  and  the  Control  Problem 

The  ultimate  objective  of  this  type  of  process  is  to  obtain  a  uniformly  distributed  metal  oxide 
deposition  of  desired  thickness  onto  a  wafer.  Currently,  the  deposition  as  function  of  the  wafer 
radius  is  roughly  exponential  with  maximum  deposition  at  the  wafer  center.  By  placing  the  wafer 
on  a  susceptor  with  larger  radius,  the  extreme  deposition  points  are  on  the  susceptor,  which 
improves  uniformity  of  the  wafer  deposition.  By  rotating  the  wafer,  the  minimum  deposition 
point  moves  from  die  edge  of  the  wafer  to  the  center,  which  also  improves  uniformity.  As  a 
result,  a  concave  deposition  profile  is  obtained. 

The  following  manipulated  inputs  are  available  to  control  deposition: 

•  flow  rate  and/or  pressure  of  the  carrier  gas; 

•  temperature  of  the  susceptor; 

•  concentration  of  the  precursor  gases. 

The  following  outputs  are  available  for  monitoring  deposition: 

•  flow  rate  measurement  of  carrier  gas  (in-situ); 

•  thermocouples  on  the  wafer  and/or  susceptor  {in-situ); 

•  deposition  thickness  measurements  (ex-situ). 

Expected  sources  of  noise  and  disturbances  are: 

•  measurement  noise  on  all  measurements; 

•  periodic  disturbance  from  wafer  rotation,  e.g.,  because  of  wafer  eccentricity; 

•  randomly  occurring  “drops”  in  flow  rate; 

•  randomly  occurring  changes  in  precursor  concentration; 

•  drifts  in  wafer  and  susceptor  temperature; 

The  control  problem  is  to  obtain  deposition  thickness  with  a  desired  deposition  uniformity  using 
the  available  control  inputs  and  measured  outputs,  in  the  face  of  the  expected  noise  and 
disturbances. 
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Figure  31.  Schematic  of  Control  Structure  for  HTS  MOCVD. 


Control  Strategy 

The  deposition  process  as  a  function  of  control  inputs  can  be  considered  as  a  static  system. 
Therefore,  there  may  be  no  need  for  a  dynamic  controller  that  tries  to  control  deposition  directly 
from  the  manipulated  inputs;  a  desired  deposition  is  obtained  by  regulating  the  manipulated 
inputs  at  a  desired  (nominal)  operating  point.  These  nominal  operating  points  have  been  identified 
in  [24].  However,  all  disturbances  are  dynanuc,  which  implies  that  dynamic  “slave”  controllers 
are  needed  locally  to  obtain  tight  regulation  at  the  desired  operating  points,  see  Figure  31.  The 
control  structure  shows  both  dynamic  a  controller  using  in-situ  sensing  as  well  as  a  run-to-run 
controller  employing  ex-situ  sensors.  The  bubblers  are  controlled  by  local  (inner-loop) 
controllers,  there  are  in-situ  temperature  sensors  measuring  substrate  temperature  corresponding 
to  substrate  (segmented)  heaters,  and  metrology  is  employed  to  measure  wafer  properties  of 
interest  (deposition  thickness  and  uniformity,  stoichiometry). 

It  is  expected  that  integral  control  and  (notch)  filtering  possibly  extended  with  some  phase 
compensation  will  give  tight  regulation.  For  this  design,  it  is  preferred  to  have  (low-complexity) 
dynamic  models  available  that  relate  measured  outputs  to  manipulated  inputs  and 
disturbances/noise. 

To  control  deposition  uniformity,  run-to-run  control  has  shown  to  be  an  excellent  means  to 
improve  current  uniformity,  see  reference  [27],  especially  the  example,  which  shows  run-to-run 
control  of  a  Rapid  Thermal  Oxidation  process  that  has  characteristics  similar  to  CVD.  The  main 
idea  is  to  adjust  the  values  of  the  nominal  operating  points  (called  recipe  variables)  from  wafer- 
to-wafer,  i.e.,  off-line  after  each  run,  see  Figure  31.  For  this  purpose,  a  static  model  is  needed 
that  relates  manipulated  inputs  to  (ex-situ)  measured  deposition  thickness.  This  model  can  be 
obtained  with  the  existing  CFD-ACE  software  by  simulating  the  deposition  process  for  varying 
values  of  one  operating  point  while  fixing  the  other  operating  points.  For  small  variations  it  is 
expected  that  only  one  variation  (fi'om  its  nominal  value)  per  operating  point  is  needed.  For 
larger  variations,  non-linear  behavior  may  be  taken  into  account  by  measuring  deposition 
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thickness  for  multiple  variations  per  operating  point.  The  (static)  map  from  operating  points 
(recipe  variables)  to  deposition  thickness  is  called  the  response  surface.  Initially,  it  is  sufficient 
to  identify  only  the  following  response  surfaces: 

•  a  nominal  response  surface  with  all  operating  points  fixed  at  their  current  nominal  values; 

•  one  response  surface  for  each  recipe  variable  with  the  operating  point  varied  from  its  nominal 
value  with  very  small  variations  (0. 1-0.5%  of  nominal  value). 

The  collection  of  response  surfaces  for  each  operating  point  is  an  w  x  n  matrix  with  m  the  number 
of  measured  deposition  thickness  points  on  the  wafer  and  n  the  number  of  recipe  variables. 
Inversion  of  this  matrix  is  the  key  to  run-to-run  control.  As  n  is  usually  much  less  than  m,  the 
maximum  obtainable  uniformity  level  is  limited.  If  this  uniformity  level  is  not  acceptable,  it  is 
possible  to  improve  it  by  increasing  the  number  of  recipe  variables,  e.g.,  by  assuming  that  there 
are  k  (independent)  temperature  control  inputs  (“rings”)  in  the  susceptor  (i.e.,  using  heater 
segmentation)  and  obtain  separate  response  surfaces  for  each  temperature  ring.  In  general,  the 
number  and  choice  of  manipulated  inputs  and  measured  outputs  can  be  a  powerful  design  variable 
to  obtain  the  desired  goal,  and  it  will  be  interesting  to  investigate  if  an  “optimal”  choice  can  be 
made. 

Run-to-Run  Control 

The  goal  of  a  manufacturing  system  is  to  produce  multiple  copies  of  the  same  product,  each 
having  attributes  within  specified  tolerances.  Product  quality  attributes  are  typically  determined 
after  the  product  is  manufactured,  since  in  most  cases  sensors  are  not  available  which  can  directly 
monitor  all  of  these  attributes  during  the  process.  Moreover,  during  the  run,  the  control  variables, 
or  recipe  variables,  are  pre-set  and  held  fixed  during  the  run.  Different  recipes  can  be  selected  to 
produce  different  products,  or  similar  products  but  with  different  attributes. 

The  run-to-run  control  problem  is  to  adjust  the  recipe  for  the  next  run  based  on  the  results  of  the 
previous  runs  such  that  the  product  quality  improves  and/or  yield  increases,  i.e.,  more  good 
product  is  produced.  This  has  proved  to  be  very  popular  in  process  control. 

In  this  section  we  describe  how  the  recipe  can  be  adjusted  from  “run-to-run”  using  a  very  simple 
algorithm  based  on  the  attributes  of  the  product  produced  in  the  previous  run  or  runs  [57].  The 
algorithm  can  be  analyzed  under  a  variety  of  assumptions  on  the  behavior  of  the  actual  process 
using  standard  control  theory  [57]. 

Let  t  =  1,2,....,  denote  the  run  number,  r,  G  R”  the  vector  of  recipe  variables  used  during  run  t,  y,  € 
R"  the  vector  of  product  quality  attributes  produced  at  the  end  of  run  t,  and  e,  the  normalized 
product  quality  error,  whose  ith  element  is  defined  as, 

^/(i)  =  Cy»(i)-  Wi)yjw/(i)  i  =  l,...,n  (1) 

where  ydesii)  is  the  ith  desired  product  quality  attribute,  and  ytoiii)  is  the  associated  error 
tolerance.  The  error,  written  in  vector  form,  becomes, 

e,=  R'^(y,-ydes) 


R  =  diag  ( ytoi(l),...,yto/(n))  (2) 

The  simplest  choice  for  run-to-run  control  is  to  correct  the  previous  recipe  by  an  amount 
proportional  to  the  current  error.  Thus,  for  run/ =1,2,...,  adjust  Ae  recipe  according  to. 
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ri=rnom  +  U, 


u,=  u,.i-Te Uo=0  (3) 

where  rnom  is  the  nominal  recipe,  m,  is  the  correction  to  the  nominal  recipe  for  run  t,  and  T  G  R™“ 
is  the  control  design  (gain)  matrix.  It  is  important  to  emphasize  that  (2)  together  with  (3) 
constitutes  the  complete  run-to-run  algorithm.  Observe  also  that  (3)  has  the  same  form  as  a 
gradient  descent  optimization  algorithm.  It  is  possible  to  choose  F  and  to  analyze  the  algorithm 
under  a  variety  of  assumptions  about  how  u,  effects  e,  [57].  It  can  be  shown  that  most  of  the 
widely  used  run-to-run  algorithms  are  in  the  form  of  (3)  for  different  choices  of  the  control 
design  (gain)  matrix  F. 

Run-to-Run  Control  for  MOCVD  Process 

In  this  section  we  will  discuss  application  of  run-to-run  control  to  the  MOCVD  process.  A 
conamon  type  of  control  problem  in  MOCVD  is  to  control  desired  stoichiometries  for  both  BaA?^ 
and  Cu/Y  by  adjusting  the  precursor  concentrations  of  Y,  Ba  and  Cu,  while  maintaining  a  desired 
growth-rate.  Hence  we  have  three  recipe  variables  available: 

'Y 
r=  Ba 
Cu 

Theoretically  we  can  have  as  many  performance  variables  as  we  wish,  but  if  this  number  is  larger 
than  the  number  of  recipe  variables,  the  computation  of  F  will  be  based  on  the  pseudo-inverse  of 
the  static  relationship  between  performance  and  recipe  variables,  and  a  singular  value 
decomposition  has  to  be  used  to  project  the  space  of  performance  variables  onto  the  space  of 
recipe  variables,  i.e.,  in  effect  only  three  performance  variables  can  be  controlled.  Consequently 
we  have  to  select  three  performance  variables.  We  choose  the  stoichiometries  of  Ba  and  Cu  at  the 
center  of  the  wafer  as  the  first  two  performance  variables,  and  the  average  growth  rate  over  all 
wafer  nodes  as  the  third: 

^Cu 
^av. 

with  Sbo  and  Scu  denoting  the  stoichiometries  of  Ba  and  Cu  at  the  center,  respectively,  and 
denoting  the  average  growth-rate  in  A/min. 

There  is  no  such  thing  as  an  optimal  stoichiometry  or  optimal  growth-rate,  rather  there  are  upper 
and/or  lower  bounds  on  the  final  stoichiometries  and  growth-rate.  For  this  process,  the  following 
desired  operating  area  has  been  identified: 

1.5  <  5b,  <2.0 
3.0<5c„<4.0 
32<R,,^40 

For  run-to-run  control  we  need  to  have  a  single  desired  operating  point,  hence  we  pick  an 
arbitrary  point  in  the  desired  area: 
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ydes 


1.8 

3.2  . 

35_ 

Furthermore,  we  have  to  compute  the  3x3  matrix  gain  F.  For  this,  we  identified  the  static 
relationship  between  perturbations  in  precursor  concentration  and  the  corresponding  perturbation 
of  the  three  performance  variables.  The  perturbations  were  taken  as  +10%  perturbations  from  the 
following  nominal  starting  points: 


r  F  1 

^  nom 

’0.000562038' 

^nom 

^^nom 

= 

0.002873710 

S'^nom  . 

_0.000327275_ 

The  corresponding  nominal  performance  variables  were  computed  as: 


nnom 

^Ba 

■  5.47229  ■ 

y  nom 

^nom 

^Cu 

= 

0.345701 

nnom 

^av 

_  33.2372  _ 

This  resulted  in  the  following  gain  matrix  : 


1.59857 

-0.370606 

1.45349 


0.675818  -0.325193 

0.611552  -0.294269 

-28.2884  -0.295681 


The  factor  is  a  number  between  0  and  1,  and  is  used  to  trade-off  speed  of  convergence  versus 
amplification  of  measurement  noise. 

We  started  the  run-to-run  iterations  with  =  However,  the  solution  diverged  rapidly  from  run 
to  run  with  these  gains.  Two  reasons  can  account  for  this.  Firstly,  the  nominal  operating  point  is 
far  away  from  the  target  operating  point.  Consequently,  the  run-to-run  gains  are  large,  because 
they  are  based  on  a  large  error.  Secondly,  the  system  is  nonlinear  while  the  run-to-run  gains  are 
based  on  a  linear  approximation.  The  combination  of  these  two  factors  make  it  impossible  to 
converge  quickly  (i.e.,  with  ^  =  1)  to  the  target  operating  point.  Therefore  we  have  to  lower  //  in 
order  to  make  the  iteration  converge.  The  downside  of  this,  however,  is  the  fact  that  (a  lot)  more 
iterations  are  required  to  converge  close  enough  to  the  target  point.  This  may  very  well  be 
undesirable,  as  more  iterations  mean  more  wafers  with  undesirable  compositions.  We  therefore 
propose  the  following  four-step  approach: 

1 .  Identify  run-to-run  gain  in  nominal  (current)  operating  point  using  (real-time)  experiments. 

2.  Simulate  a  run-to-run  iteration  with  an  approximate  (linear)  system  model  using  a  low  value 
of  //  and  use  as  many  iterations  so  as  to  converge  with  high  accuracy  to  the  desired  operating 
point. 

3.  Use  the  optimal  simulated  input  values  as  a  starting  point  for  a  new  experiment.  Determine 
how  much  the  measured  performance  differs  from  the  optimal  simulated  performance. 

4.  If  necessary,  identify  new  run-to-run  gains  in  the  new  operating  point  -  which  should  be 
much  closer  to  the  target  operating  point  -  and  repeat  the  run-to-run  experiment. 
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We  applied  this  four-step  approach  to  the  MOCVD  process.  With  the  identified  gain  f  we  started 
step  2,  the  run-to-run  simulation.  Figures  32  and  33  show  thirty  run-to-run  iterations  in  a  ternary 
diagram  for  two  different  values  of  n,  namely  fx  =  0.2  and  =  0.3,  respectively.  The  iteration 
with  n  =  0.3  is  almost  unstable,  and  larger  values  of  n  indeed  make  the  iteration  diverge. 
Iterations  with  =  0.1  (not  shown)  look  similar  to  iterations  with  [x  =  0.2,  except  that  now 
approximately  50  iterations  are  needed  to  reach  the  target  point.  Clearly,  a  run-to-run  iteration 
with  fx  =  0.2  converges  to  the  target  operating  point  with  any  desired  accuracy.  Also,  the  average 
growth-rate  converged  to  its  desired  value  of  35.  After  thirty  run-to-run  iterations,  the 
performance  variables  converged  to: 


30 

3’//=0.2  = 

'1.95305' 

3.19925 

30 

.  y;/=0.3  = 

'1.84445' 

2.93252 

35.5246 

35.1611 

We  recomputed  the  simulation  with  fx  =  0.2,  now  performing  200  iterations  so  that  the  converged 
values  reached  the  target  value  within  an  accuracy  of  le-5.  The  resulting  optimal  input  values  are 
given  by: 


,.200  _ 
0/=0.2  - 


0.000931721 

0.001566990 

0.005022070 


With  these  values  we  performed  step  three  of  the  four-step  approach.  We  recomputed  the 
nonlinear  CFD-ACE  simulation  to  check  how  much  this  solution  would  differ  from  the 
approximate  solution.  The  resulting  performance  variables  were  given  by: 


3'  = 


1.80745' 

3.16043 

35.0786 


which  is  very  close  to  the  desired  values  of  1.8,  3.2  and  35,  respectively.  This  shows  the  validity 
of  the  four-step  approach.  In  effect,  only  five  ‘expensive’  experiments  had  to  be  conducted  in 
order  to  converge  with  high  accuracy  to  the  target  operating  point. 

Figure  34  checks  how  controlling  the  stoichiometries  at  the  center  point  and  controlling  the 
average  growth-rate  affect  the  stoichiometries  and  the  growth-rate  at  the  other  wafer  nodes, 
respectively.  Apparently,  the  values  of  the  stoichiometries  and  the  growth-rate  at  each  wafer 
node  lie  within  the  desired  operating  range.  Given  the  fact  that  the  shape  of  the  curve  is  very 
deterministic,  one  could  similarly  control  one  of  the  extreme  points  at  the  wafer  edge  or  one  of 
the  points  between  the  edge  and  the  center,  and  check  the  other  points  afterwards.  A  motivation 
for  choosing  another  point  than  the  center  point,  might  be  that  the  other  points  £ire  closer  to  the 
target  point,  hence  convergence  of  the  iteration  will  be  faster.  If  one  or  more  wafer  nodes  fall 
outside  the  desired  operating  range,  one  can  change  the  target  operating  point  and  shift  it  to  a 
position  inside  the  desired  area  where  all  wafer  nodes  meet  the  specification,  or  one  can  control 
another  wafer  node.  If  desired,  one  can  control  different  wafer  nodes  for  Ba  and  Cu 
stoichiometries. 
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Ternary  Diagram 


Cu 


Figure  32.  Ternary  diagram  showing  30  run-to-run  iterations  from  the  nominal  starting  point  to 
the  target  operating  point  with  //  =  0.2.  The  starting  point  and  target  point  are  shown  as  a  thick 
green  and  red  bullet,  respectively,  while  the  iterations  are  shown  as  small  triangles.  Also  shown 
in  red  is  the  desired  operating  area. 


Ternary  Diagram 
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Figure  33.  Ternary  diagram  showing  30  run-to-run  iterations  from  the  nominal  starting  point  to 
the  target  operating  point  with  /r  =  0.3.  The  starting  point  and  target  point  are  shown  as  a  thick 
green  and  red  bullet,  respectively,  while  the  iterations  are  shown  as  small  triangles.  Also  shown  in 
red  is  the  desired  operating  area. 
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Figure  34.  Final  performance  for  all  wafer  nodes  along  the  wafer  diameter,  (a)  Ba  stoichiometry; 
(b)  Cu  stoichiometry;  (c)  growth-rate. 

The  implication  of  this  result  is  very  significant.  It  suggests  that  it  is  possible  to  speed  up  the  run- 
to-run  control  by  performing  the  run-to-run  iterations  in  the  simulation  (rather  than  experiment) 
and  only  apply  control  to  the  actual  chamber  once  the  stoichiometry  is  within  the  desired  region. 
The  run-to-run  control  iterations  via  the  virtual  chamber  avoid  the  cost  associated  with  the  use  of 
actual  wafers  and  increase  throughput  significantly  as  well  as  improving  reproducibility.  Even  if 
the  run-to-run  control  cannot  provide  the  desired  stoichiometry  on  the  first  try,  it  will  provide  the 
ability  to  get  much  closer  to  the  desired  stoichiometry  region  than  would  otherwise  be  possible. 
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6.  Conclusions 

We  have  developed  a  systematic  methodology  to  look  at  the  stracture  and  reactivity  of  the  metal 
P-diketonates  precursors,  used  in  the  chemical  vapor  deposition  of  high  temperature 
superconductors,  using  first  principle  quantum  chemistry  methods.  The  large  size  of  the  precursor 
molecules  has  so  far  limited  studies  of  this  nature.  Very  little  information  is  available  on  the 
chemistry  of  these  compounds.  Semi-empirical  methods  were  found  inadequate  in  describing  the 
structures  and  traditional  ab-initio  methods  were  deemed  to  be  very  expensive,  given  the  size  of 
the  molecules.  We  have  tested  and  validated  different  density  functionals  for  describing  these 
molecules.  Based  on  our  validations,  an  effective  core  potential  was  chosen  to  describe  the  metal 
(Cu  and  Y).  We  have  investigated  the  decomposition  mechanism  for  Cu(acac)2  and  Cu(dpm)2, 
two  of  the  most  common  precursors  for  copper  in  HTS  films.  A  kinetic  mechanism  for 
decomposition  was  formulated  and  the  initial  break-up  of  the  Cu-0  bond  was  determined  to  be 
the  rate-limiting  step.  A  similar  investigation  was  also  carried  out  for  Y(acac)2,  which  serves  as  a 
model  for  Y(dpm)3.  A  preliminary  mechanism  has  been  put  together  based  on  the  work  done  so 
far. 

Results  have  been  obtained  for  reactor-scale  physical  modeling  of  MOCVD  of  YBCO.  The 
model  includes  chemical  kinetics  and  species  transport  to  the  wafer  surface.  The  YBCO 
deposition  rates  obtained  from  these  simulations  are  comparable  to  those  published  in  the 
literature.  Sensitivity  calculations  using  this  model  were  used  for  run-to-run  controller  design. 

A  general  control  structure  was  developed  for  MOCVD  reactor  control.  An  innovative  run-to-run 
control  method  was  developed  which  enables  efficient  stoichiometry  control.  A  standard  run  is 
performed  and  performance  measures  are  determined  using  ex-situ  metrology.  The  run-to-run 
control  iterations  are  then  run  on  the  virtual  reactor  model  to  suggest  control  perturbations  for  the 
next  actual  run  to  reach  the  desired  stoichiometry  region. 
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7.  Future  Work 

Further  work  needs  to  be  done  to  get  a  more  detailed  picture  of  the  decomposition  pathways.  In 
particular,  the  effect  of  the  oxidant  gas  such  as  O2  and  N2O,  if  any,  needs  to  be  explored  further. 
In  this  study,  we  have  concentrated  our  attention  on  copper  precursors  due  to  the  computational 
efforts  involved  in  study  of  larger  molecules  such  as  yttrium  and  barium.  A  similar  detailed 
treatment  needs  to  be  carried  out  for  yttrium  precursors  and  eventually  extended  to  barium 
precursors.  Yttrium  precursors  are  known  to  coordinate  a  water  molecule  to  give  Y(dpm)3.H20 
molecule  [1].  The  effect  of  this  coordination  on  the  decomposition  pathways  needs  to  be 
explored  further. 

The  detailed  models  need  to  be  reduced  for  reactor  modeling,  simulation,  and  control  design.  We 
will  continue  investigating  the  application  of  POD  and  other  model-order  reduction  methods  to 
the  full  chemical-transport  model. 

The  refined  kinetic  mechanism  obtained  from  the  yttrium  and  barium  calculations  need  to  be 
incorporated  into  the  reactor  scale  model.  The  predictions  of  the  reactor  scale  model  need  to  be 
validated  against  the  experimental  results  such  as  those  planned  as  part  of  the  Caltech  VIP  effort. 

The  procedures  used  in  this  work  could  be  expanded  to  ferroelectrics,  such  as  BaSrTiOs,  since 
similar  MOCVD  precursors  are  used  [58].  The  chemistry  of  MOCVD  of  ferroelectrics  is  not  well 
understood,  yet  there  is  a  considerable  interest  in  using  these  materials  in  ULSI IC  production.  A 
detailed  study  of  the  basic  reaction  pathways  and  rates,  similar  to  the  one  performed  for  HTS, 
would  provide  the  necessary  data  for  simulation  of  MOCVD  of  ferroelectrics. 
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8.  Appendix  L  Data  for  the  Quantum  Chemistry  Validation 
Study 


Cu 


Excitation  energies  Eexec  and  ionization  potentials  BP  for  Cu  2S  with  different  basis  sets  and 
methods  [eV] _ 


Eexec 

I.P. 

ECPl 

CAS 

-2.88 

B3LYP 

8.00 

7.94 

CCSD(T)  (3d/4s)* 

1.69 

6.91 

ECP2 

CAS 

-3.00 

B3LYP 

6.77 

8.68 

CCSD(T)  (3d/4)* 

1.32 

6.88 

CCSD(T) 

1.82 

6.85 

CCSD(T)  ^ 

1.65 

7.13 

ECP-S[28] 

CAS 

-3.12 

B3LYP 

7.01 

8.66 

CCSD(T)  (3d/4s)* 

1.47 

7.07 

CCSD(T) 

1.76 

7.04 

Experimental 

[46] 

1.49 

7.73 

“  includes  only  correlation  between  the  3d  and  4s/p  shell 
**  additional  d  and/function 


Ionization  potentials  for  Cu  2S  with  different  basis  sets  and  methods  [eV] 


BLYP 

B3LYP 

LSDA 

6-3 IIG 

4.15 

6.86 

8.06 

6-31 1+G 

8.18 

8.02 

8.82 

6-3 IIG** 

7.78 

7.48 

8.06 

6-31 1+G** 

8.18 

8.02 

8.82 

ECPl 

7.67 

7.67 

8.27 

ECP2 

9.33 

7.82 

8.53 

CuH 

Method 

Basis  Set 

Te 

CO 

De[eV] 

Cl 

aU-electron 

1.59 

1584 

2.36  [47] 

CI(cv) 

1.50 

1740 

2.49 

CI(cc) 

1.51 

1815 

2.42 

B31YP 

ECPl 

1.51 

1868 

3.29 

ECP2 

1.50 

1837 

2.73 

631 1+G** 

1.48 

1879 

2.69 

Exp. 

1.46 

1941 

2.85 

CuMe,  CuEt 


Dissociation  energies  [De]  for  CuMe  and  CuEt  [kcal/mol] 
CuMe 

experimental:  about  60  kcal/mol 

other  theoretical  work:  Bauschlicher  [48]  45.8 


Ziegler  [49]  56.8 
Boehme  [31]  58.0 


Method  Basis  Set 

CuMe  CuEt 

B3LYP  6-3 IIG 

78.82 

6-311+G 

45.25 

6-3 IIG** 

65.78 

6-311+G** 

52.37  45.98 

AhlrichsTZV 

ECPl 

59.09  53.80 

ECP2 

53.66  48.01 

LSDA  6-3 IIG 

6-311+G 

6-3 IIG** 

6-311+G** 

78.17  70.89 

ECPl 

90.59  84.64 

ECP2 

81.04  74.82 

BLYP  6-3 IIG 

6-311+G 

6-3 IIG** 

6-311+G** 

56.96  50.85 

ECPl 

69.09  63.78 

ECP2 

25.71  20.49 

BSSEforCuM 

e  using  B3LYP  [kcal/mol] 

Basis  Set 

BSSE  corrected  De  uncorrected  De 

6-311+G** 

ECPl 

ECP2 

-5.56  57.93  52.37 

7.95  51.14  59.09 

-2.89  56.55  53.66 

Comparison  of  calculated  and  experimental  CuMe  frequencies 


B3LYP 

calculated  frequencies 

6-311+G**  ECPl 

ECP2 

Experimental 

vibrations 

theo/exp 

6-311+G** 

ECPl 

ECP2 

Cu-Me  str 

522 

510 

533 

350 

1.49 

1.46 

1.52 

Cu-Me  tilt 

674 

735 

692 

424 

1.59 

1.73 

1.63 

CH  sy-b 

1122.7 

1293.8 

1152 

CH  as-b 

1445.2 

1600.7 

1469 

1336 

1.08 

1.20 

1.10 

CH  as-str 

3126.3 

3447.4 

3158 

2929 

1.07 

1.18 

1.08 

CH  sy-str 

3030.6 

3269.6 

3036 

2880 

1.05 

1.14 

1.05 

BLYP 

6311pss 

MB 

DZ 

Cu-Me  str 

516 

507 

530 

1.47 

1.45 

1.51 

Cu-Me  tilt 

671 

722 

683 

1.58 

1.70 

1.61 

54 


CH  sy-b 

1092 

1248 

1120 

CH  as-b 

1406 

1547 

1426 

1.05 

1.16 

1.07 

CH  as-str 

3054 

3334 

3082 

1.04 

1.14 

1.05 

CH  sy-str 

2959 

3159 

2961 

1.03 

1.10 

1.03 

LSDA 

6311pss 

MB 

DZ 

Cu-Me  str 

590 

586 

603 

1.69 

1.67 

1.72 

Cu-Me  tilt 

648 

737 

685 

1.53 

1.74 

1.62 

CH  sy-b 

1060 

1260 

1115 

CH  as-b 

1358 

1535 

1397 

1.02 

1.15 

1.05 

CH  as-str 

3083 

3407 

3116 

1.05 

1.16 

1.06 

CH  sy-str 

2971 

3228 

2981 

1.03 

1.12 

1.04 

Calculated  harmonic  fundamentals  for  CuEt 


Calculated  harmonic  fundamentals  for  CuEt  (relative  to  6-3 I  IG) 


B3LYP 


6311+G 


0.98 


0.86 


0.89 


0.93 


6311G** 


0.99 


1.00 


0.99 


0.95 


0.97 


0.97 


1.00 


0.97 


0.98 


0.97 


631 1+G** 

ECPI 

ECPn 

0.97 

0.96 

0.97 

0.86 

0.86 

0.83 

0.88 

0.86 

0.90 

0.90 

0.96 

0.92 

0.97 

1.06 

0.98 

0.96 

1.05 

0.98 

0.99 

1.09 

1.00 

6.93 

1.03 

0.95 

0.97 

1.05 

0.98 

0.97 

1.07 

0.98 
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1.00 

0.98 

0.98 

1.08 

0.99 

1.00 

0.97 

0.97 

1.09 

0.99 

CH  b  asym 

1.00 

0.97 

1.09 

0.99 

CHS  str  einzel 

1.00 

ieshhh 

1.00 

1.10 

1.01 

1.00 

1.01  • 

1.01 

1.10 

1.01 

CH2  str  asym 

1.00 

1.00 

1.00 

1.12 

1.01 

CHS  str  symm 

1.00 

1.01 

1.01 

I.IS 

1.02 

CHS  str  asy 

1.00 

1.00 

1.00 

I.IS 

1.02 

CuO 


Dissociation  energies  De  [kcal/mol] 


Method 

Basis  Set 

DeCuO 

LSDA 

6-SllG 

I1S.4 

6-Sll+G 

II8.1 

6-SllG** 

92.5 

6-Sll+G** 

97.1 

ECPI 

84.0 

ECPn 

91.1 

BLYP 

6-SllG 

170.4 

6-Sll+G 

17S.S 

6-SllG** 

70.0 

6-Sll+G** 

72.6 

ECPI 

66.7 

ECPn 

70.S 

BSLYP 

6-SllG 

75.S 

6-Sll+G 

78.8 

6-SllG** 

59.1 

6-Sll+G** 

61.7 

ECPI 

45.0 

ECPn 

58.0 

Experimental[50] 

62.7  V  S.O 
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9.  Appendix  II:  Publications 

Journal  Papers: 


1.  “A  Density  Functional  Theory  Study  of  (Hfac)Cu-L  Compounds  of  interest  for  the  Chemical 
Vapor  Deposition  of  Copper,”  Carlo  Cavallotti,  Vijay  Gupta,  Cornelia  Sieber  and  Klavs  F. 
Jensen,  to  be  submitted  to  Journal  of  Physical  Chemistry. 

2.  “Computational  Chemistry  Investigations  of  Precursor  Decomposition  Pathways  in  MOCVD 
of  High  Temperature  Superconductors,”  Vijay  Gupta  and  Klavs  F.  Jensen,  to  be  submitted  to 
Journal  of  Crystal  Growth. 

Conference  Presentation: 

1.  S.  Ghosal,  J.  L.  Ebert,  D.  de  Roover,  and  A.  Emami-Naeini,  “Model-based  Control  of 

MOCVD  Rate,  Uniformity  and  Stoichiometry,”  Presented  at  the  Third  Symposium  of  Process 
Control,  Diagnostics,  and  Modeling  in  Semiconductor  Manufacturing,  795'*  Meeting  of  the 
Electrochemical  Society,  Seattle,  May  2—6, 1999. 
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